Stochastic Modeling of Particle Transport in Confined Geometries: Problems and Peculiarities

Abstract

The equivalence between parabolic transport equations for solute concentrations and stochastic dynamics for solute particle motion represents one of the most fertile correspondences in statistical physics originating from the work by Einstein on Brownian motion. In this article, we analyze the problems and the peculiarities of the stochastic equations of motion in microfluidic confined systems. The presence of solid boundaries leads to tensorial hydrodynamic coefficients (hydrodynamic resistance matrix) that depend also on the particle position. Singularity issues, originating from the non-integrable divergence of the entries of the resistance matrix near a solid no-slip boundary, determine some mass-transport paradoxes whenever surface phenomena, such as surface chemical reactions at the walls, are considered. These problems can be overcome by considering the occurrence of non vanishing slippage. Added-mass effects and the influence of fluid inertia in confined geometries are also briefly addressed.

1. Introduction

The main legacy of the Einsteinian theory of Brownian motion to modern physics lies in the confirmation of the atomistic nature of matter and in the equivalence between random molecular motion at the microscopic level and the macroscopic phenomenon of diffusion [1,2].

Thanks also to the contributions by Langevin, Smoluchowski and many others [3,4,5], it is now common knowledge in applied transport theory [6] that any transport equation, expressed in the form of an advection-diffusion equation for the concentration field $c(\mathbf{x},t)$ of some solute diffusing in a fluid phase,

$$\frac{\partial c(\mathbf{x},t)}{\partial t}=-\nabla ·\left[{\mathbf{u}}^{\left(p\right)}\left(\mathbf{x}\right)c(\mathbf{x},t)\right]+D{\nabla }^{2}c(\mathbf{x},t)$$

(1)

where ${\mathbf{u}}^{\left(p\right)}\left(\mathbf{x}\right)$ is the macroscopic velocity field experienced by the solute (coinciding in most of the applications with the fluid phase velocity $\mathbf{u}\left(\mathbf{x}\right)$), and D its isotropic diffusivity, can be represented in terms of the microscopic motion of the solute particles by means of a Langevin equation of the form

$$d\mathbf{x}\left(t\right)={\mathbf{u}}^{\left(p\right)}\left(\mathbf{x}\left(t\right)\right)dt+\sqrt{2D}d\mathbf{w}\left(t\right)$$

(2)

where $d\mathbf{w}\left(t\right)=(d{w}_1\left(t\right),d{w}_2\left(t\right),d{w}_3\left(t\right))$ is the increment of a three-dimensional Wiener process in the time interval $(t,t+dt)$. This equivalence is also computationally important as it enables to solve parabolic transport equations in complex geometries by means of stochastic simulations of the Langevin Equation (2) [7,8,9,10].

The description of physical processes in terms of stochastic Langevin equations, both in equilibrium and in out-of equilibrium conditions, has become one of the strongest and more fruitful research lines in modern statistical physics, providing useful insights in all the fields of physical investigation, including quantum and particle physics [11,12,13].

The Langevin Equation (2), expressed exclusively with respect to the particle position $\mathbf{x}\left(t\right)$, can be derived from a stochastic dynamics (stochastic Newton equation) of the form

$$\begin{array}{ccc}\hfill \frac{d\mathbf{x}}{dt}& =& \mathbf{v}\hfill \\ \hfill m\frac{d\mathbf{v}}{dt}& =& -\eta \left[\mathbf{v}-{\mathbf{u}}^{\left(p\right)}\left(\mathbf{x}\right)\right]+{\mathbf{F}}_{\mathrm{stocha}}\left(t\right)\hfill \end{array}$$

(3)

where ${\mathbf{F}}_{\mathrm{stocha}}\left(t\right)dt=\alpha d\mathbf{w}\left(t\right)$, $\alpha =\sqrt{2k_{b}T\eta }$, m is the particle mass, $\eta$ the friction factor, T the absolute temperature and $k_B$ the Boltzmann constant, in the limit for $m/\eta \ll 1$, i.e., in the case particle inertia could be neglected. This simplification is not only important in transport theory, but also in the thermodynamics of fluctuations as it implies the simplest form of fluctuation-dissipation relation [14,15],

$$D\eta =k_{B}T$$

(4)

that, in the case of spherical particles for which $\eta =6\pi \mu R_p$ ($\mu$ is the fluid viscosity, $R_p$ the particle radius) is referred to as the Stokes-Einstein relation. Equation (4) connects a transport parameter, related to the intersity of fluctuations (the diffusivity D) to a hydrodynamic quantity, related to dissipation (the friction factor $\eta$) at constant temperature T.

The rise of microfludics as an autonomous hydrodynamic discipline has informed recent research in the theory of fluid motion, with relevant practical applications in analytical chemistry, molecular biology and biomedicine, production of high-quality pharmaceuticals, etc. [16,17,18,19]

The analysis of flow systems at microscale in devices (microchannels, microreactors) of almost comparable size than that of solute particles, necessarily implies a more detailed understanding of all the physical processes at small lengthscales, and thus a unitary and integrated use of theoretical tools deriving from practically all the branches of physics coupled to hydrodynamic problems [20]: optical methods for trapping, particle manipulation and detection [21,22] are of common use in microfluidics, and electrosmotic [23,24,25] and acoustic effects [26,27], either for fluid mixing or for separation purposes can be conveniently integrated within a microfluidic device. On the other hand, microscale hydrodynamics represents the natural realm where thermal and hydrodynamic fluctuations play a leading role. This admits also interesting practical implications related to the use of micrometric particles for detecting and measuring superficial phenomena and surface interactions (Brownian probes) [28,29,30]. For 1 $\mathsf{\mu }$m particles, quantum effects and forces become measurable, and this opens the way for testing quantum theories, such as the occurrence of the Casimir effect, using “macroscopic” microfluidic systems [31,32], paving the way for a better understanding of the quantum-to-classical transition.

In most of the microfluidic applications, linearized hydrodynamic equations can be used with reasonable accuracy [33,34]. For timescales much larger than the dissipation time $t_{\mathrm{diss}}=m/\eta$, the instantaneous Stokes equations

$$\left\{\begin{array}{c}\mu {\nabla }^2\mathbf{u}\left(\mathbf{x}\right)-\nabla p\left(\mathbf{x}\right)=-\nabla ·\mathsf{\tau }\left(\mathbf{x}\right)=0\hfill \\ \nabla ·\mathbf{u}\left(\mathbf{x}\right)=0\hfill \end{array}\right.$$

(5)

where $(\mathbf{u},p,\mathsf{\tau })$ are the velocity, the pressure and the stress tensor of the fluid respectively, provides a good approximation for modeling fluid-particle interactions, and the no-slip boundary conditions at the particle surface ($S_p$)

$$\mathbf{u}\left(\mathbf{x}\right)=\mathbf{v}+(\mathbf{x}-{\mathbf{x}}_0)\times \mathsf{\omega }\mathbf{x}\in S_p$$

(6)

where $(\mathbf{v},\mathsf{\omega })$ is the velocity and the angular velocity of the particle, and ${\mathbf{x}}_0$ a reference point, say the center of mass of the particle, describe with sufficient accuracy fluid-particle interactions.

Consequently, all the theory of Stokesian hydrodynamics [33,34,35,36], developed throughout more than one and half century, finds in microfluidics a vast field of application. Exceptions to the applicability of the Stokesian hydrodynamics in microfluidics occur for specific inertial effects such as the Segre-Silberberg effect [37,38,39], and the use of “high-Reynolds number flows” (with a Reynolds number order of 100–1000) in microchannels (T- and Y-junctions) for nanoparticle production [40,41,42].

The linearity of the Stokes equations and of the boundary conditions at the fluid-particle interface yields a the linear relation between the velocity of the particle and the hydrodynamic force acting on it, that in the case of a sphere moving in an unbounded fluid consist in the well known Stokes’ laws, $\mathbf{F}=-{\eta }_{\infty }\mathbf{v}$ and $\mathbf{T}=-{\eta }_{\infty }^{\omega }\mathsf{\omega }$, ($\mathbf{F}$ is the force and $\mathbf{T}$ the torque), where ${\eta }_{\infty }=6\pi \mu R_p$ and ${\eta }_{\infty }^{\omega }=8\pi \mu R_p^3$ are scalar coefficients due to the isotropy of the problem and the uniformity of the fluid.

On the other hand, whenever the fluid is confined into a bounded device, corresponding to the hydrodynamic domain bounded by walls, both the symmetries (isotropy and uniformity) are generically broken. Therefore, the hydrodynamic resistance law attains a more general position-dependent and tensorial character

$$\left(\begin{array}{c}\mathbf{F}\\ \mathbf{T}\end{array}\right)=-\mathbf{H}\left(\mathbf{x}\right)\left(\begin{array}{c}\mathbf{v}\\ \mathsf{\omega }\end{array}\right)$$

(7)

where $\mathbf{H}\left(\mathbf{x}\right)$ is the $6\times 6$ overall resistance matrix of the hydrodynamic interactions, possessing a block structure

$$\mathbf{H}=\left(\begin{array}{ccc}\mathsf{\eta }\left(\mathbf{x}\right)& & {\mathbf{C}}^{\left(1\right)}\left(\mathbf{x}\right)\\ {\mathbf{C}}^{\left(2\right)}\left(\mathbf{x}\right)& & {\mathsf{\eta }}^{\omega }\left(\mathbf{x}\right)\end{array}\right)$$

(8)

where (i) $\mathsf{\eta }$, and ${\mathsf{\eta }}^{\omega }$ are the translational and rotational friction matrices, respectively, (ii) ${\mathbf{C}}^{\left(1\right)}$, and ${\mathbf{C}}^{\left(2\right)}$ are the roto-translational coupling matrices. By the Lorentz’s reciprocal theorem it is possible to prove [33] that $\mathsf{\eta }$ and ${\mathsf{\eta }}^{\omega }$ are symmetric matrices, while ${\mathbf{C}}^{\left(1\right)}$ and ${\mathbf{C}}^{\left(2\right)}$ satisfy the property ${\mathbf{C}}^{\left(2\right)}={\left[{\mathbf{C}}^{\left(1\right)}\right]}^T$, where the superscript “T” indicates the transpose. Correspondingly, the overall resistance matrix $\mathbf{H}\left(\mathbf{x}\right)$ is symmetric and positive definite.

This automatically implies, due to the fluctuation-dissipation relation (4), a tensorial and position-dependent diffusivity. The latter two properties have deep and non trivial implications whenever a stochastic equation of motion, in the form of Equation (3) is considered, due to the highly singular nature of the Wiener description of the thermal/hydrodynamic fluctuations. For the sake of completeness, it should be also mentioned that a position dependent effective diffusivity arises in modeling solute transport in microchannels with undulated walls, in the case the transport problem is referred exclusively to the channel axial coordinate [43,44]. This is referred to as the Fick-Jacobs approximation and it is essentially a geometrical effect within an approximate transport model unrelated to any hydrodynamic interactions.

The scope of this article is to analyze in detail the problems and the peculiarities of particle transport in microfluidic systems, originating from the confined nature of the flow, starting from the stochastic description of the microscopic particle motion, in a way that may also be useful for researchers in microfluidics that are not fully familiar with stochastic differential equations. While stochastic modeling of particle transport involving constant effective diffusivities is widely used in the analysis of microfluidic devices, the inclusion of hydrodynamic effects, deriving from fluid confinement, represents a completely new and unexplored field of theoretical and numerical investigation. This article attempts to fill this gap, providing a methodological bridge between hydrodynamic theory and the stochastic formulation of transport phenomena in confined geometries, addressing also in a clear and critical way the complexities and the difficulties of this approach. The article addresses also novel original derivations as regards specific topics, such as the use of the overdamped approximation in non-equilibrium conditions and the effect of slippage as regards transport phenomena involving surface effects (such as surface chemical reactions).

In order to illustrate the physical concepts and the resulting analytical approaches, the hydrodynamics of a spherical particle near an infinite planar solid surface is explicitly considered, as a prototypical and paradigmatic case study of the hydrodynamic problems occurring in microfluidic channels. Throughout this article, rigid particles and Newtonian fluids are considered.

The article is organized as follows. Section 2 addresses the formulation of fluctuation dissipation relations in confined systems, and the computational problems associated with it. Section 3 analyzes the reduction of the equation of motion in the form of a Langevin Equation (3), discussing also the case of non-equilibrium thermal conditions (thermophoresis). Section 4 discusses the problems arising from the non-integrable singularity of the friction factor near a solid no-slip wall, as regards diffusional problems in the presence of superficial phenomena (surface chemical reactions). A way for overcoming these problems lies in relaxing the no-slip assumption, as discussed in Section 5, where we show the necessity of considering slippage effects both at the particle external surface and at the walls of the fluid domain. Finally, Section 6 discusses in a succinct way the role of the fluid inertial effects in the formulation of the stochastic equations for particle motion.

2. Fluctuation-Dissipation Relations in Confined Geometries

Consider the motion of a rigid spherical particle, of mass m and momentum of inertia I, in a microfluidic device. Let $\mathbf{v}$ and $\mathsf{\omega }$ be its velocity and angular velocity, respectively. Assume that the particle is dispersed in a fluid flow, characterized by the fluid velocity $\mathbf{u}\left(\mathbf{x}\right)$ and that the particle is also subjected to the action of an external potential field $U\left(\mathbf{x}\right)$. The particle equations of motion read

$$\begin{array}{ccc}\hfill md\mathbf{v}& =& ({\mathbf{F}}_{\mathrm{flow}}+{\mathbf{F}}_{\mathrm{hydro}})dt-\nabla Udt+{\mathbf{F}}_{\mathrm{stocha}}dt\hfill \\ \hfill Id\mathsf{\omega }& =& ({\mathbf{T}}_{\mathrm{flow}}+{\mathbf{T}}_{\mathrm{hydro}})dt+{\mathbf{T}}_{\mathrm{stocha}}dt\hfill \end{array}$$

(9)

where ${\mathbf{F}}_{\mathrm{flow}}$ and ${\mathbf{T}}_{\mathrm{flow}}$ are the force and the torque deriving from the action of the external flow field $\mathbf{u}\left(\mathbf{x}\right)$ (assuming that the particle velocity and angular velocity are vanishing), ${\mathbf{F}}_{\mathrm{hydro}}$ and ${\mathbf{T}}_{\mathrm{hydro}}$ are the force and the torque due to the hydrodynamic interactions (described by means of Equation (7) in the case the external flow velocity is vanishing and the particle possesses a velocity $\mathbf{v}$ and an angular velocity $\mathsf{\omega }$), and ${\mathbf{F}}_{\mathrm{stocha}}$, ${\mathbf{T}}_{\mathrm{stocha}}$ represent the the stochastic force and the torque deriving from to thermal fluctuations. This decomposition is made possible because the hydrodynamic equations for the fluid are assumed to be linear. The basic problem in the statistical physics of microparticle motion resides in the determination of the contributions of the thermal perturbations ${\mathbf{F}}_{\mathrm{stocha}}$ and ${\mathbf{T}}_{\mathrm{stocha}}$, since all the other terms in Equation (9) stem from a classical hydrodynamic analysis.

Following the original approach due to Einstein and Langevin [15], in the case the fluid is described by means of an instantaneous response (Stoke’s regime), it is natural to represent ${\mathbf{F}}_{\mathrm{stocha}}$ and ${\mathbf{T}}_{\mathrm{stocha}}$ in the form of a linear superposition of vector-valued Wiener processes, i.e., as

$$\begin{array}{c}\hfill {\mathbf{F}}_{\mathrm{stocha}}\left(\mathbf{x}\right)dt=\mathsf{\alpha }\left(\mathbf{x}\right)d\mathbf{w}\left(t\right)+\mathsf{\gamma }\left(\mathbf{x}\right)d{\mathbf{w}}^{\omega }\left(t\right)\\ \hfill {\mathbf{T}}_{\mathrm{stocha}}\left(\mathbf{x}\right)dt=\mathsf{\epsilon }\left(\mathbf{x}\right)d\mathbf{w}\left(t\right)+\mathsf{\beta }\left(\mathbf{x}\right)d{\mathbf{w}}^{\omega }\left(t\right)\end{array}$$

(10)

where $d\mathbf{w}\left(t\right)=(d{w}_1\left(t\right),d_2\left(t\right),d{w}_3\left(t\right))$ and $d{\mathbf{w}}^{\omega }\left(t\right)=(d{w}_1^{\omega }\left(t\right),d{w}_2^{\omega }\left(t\right),d{w}_3^{\omega }\left(t\right))$ are the increments in the time interval $(t,t+dt)$ of two mutually independent vector-valued Wiener processes. This observation is a consequence of the fact that Wiener processes are also memoryless, in the meaning that if one defines $\mathsf{\xi }\left(t\right)=d\mathbf{w}\left(t\right)/dt=({\xi }_1\left(t\right),{\xi }_2\left(t\right),{\xi }_3\left(t\right))$, interpreted in a distributional meaning, than the correlation function $\langle {\xi }_i(t_0,+t){\xi }_j\left(t_0\right)\rangle =\delta \left(t\right){\delta }_{i,j}$ is impulsive (here $\langle ·\rangle$ indicates indifferently either ensemble or temporal averages, and $t_0>0$ is any time instant) [45].

Henceforth, in order to simplify the notation, the explicit dependence of the matrices $\mathsf{\alpha },\mathsf{\gamma },\mathsf{\epsilon },\mathsf{\beta }$ on the position $\mathbf{x}$ will be omitted. While the determination of ${\mathbf{F}}_{\mathrm{flow}}$, ${\mathbf{F}}_{\mathrm{hydro}}$ and ${\mathbf{T}}_{\mathrm{flow}}$, ${\mathbf{T}}_{\mathrm{hydro}}$ follows for the simple application of the Stokesian hydrodynamics, the estimate of the matrices entering Equation (10) and defining the thermal perturbations requires a statistical physical ansatz that, at constant temperature T, the thermal fluctuations described by Equation (10) would provide the known result of equilibrium statistical physics [14]. This is the essence of the fluctuation-dissipation ansatz, and for this reason, owing to linearity, it is sufficient to consider the statistical properties for the particle dynamics in the absence of external forcings, i.e., for ${\mathbf{F}}_{\mathrm{flow}}=\nabla U=0$, ${\mathbf{T}}_{\mathrm{flow}}=0$. Under these conditions, by substituting Equations (7), (8) and (10) into the equations of motion (9) one obtains,

$$\begin{array}{ccc}\hfill d\mathbf{v}\left(t\right)& =& -\frac{\eta }{m}\mathbf{v}\left(t\right)dt-\frac{{\mathbf{C}}^{\left(1\right)}}{m}\mathsf{\omega }\left(t\right)dt+\frac{\alpha }{m}d\mathbf{w}\left(t\right)+\frac{\gamma }{m}d{\mathbf{w}}^{\omega }\left(t\right)\hfill \\ \hfill d\mathsf{\omega }\left(t\right)& =& -\frac{{\mathbf{C}}^{\left(2\right)}}{I}\mathbf{v}\left(t\right)dt-\frac{{\eta }^{\omega }}{I}\mathsf{\omega }\left(t\right)dt+\frac{\epsilon }{I}d\mathbf{w}\left(t\right)+\frac{\beta }{I}d{\mathbf{w}}^{\omega }\left(t\right)\hfill \end{array}$$

(11)

where also for $\eta ,{\mathbf{C}}^{\left(1\right)},{\mathbf{C}}^{\left(2\right)},{\mathsf{\eta }}^{\omega }$ the explicit dependence on the position has been omitted.

The matrices $\mathsf{\alpha }$, $\mathsf{\beta }$, $\mathsf{\gamma }$, $\mathsf{\epsilon }$ are the stochastic amplitude matrices, and the final goal of fluctuation-dissipation analysis is their determination from physical principles enforcing equilibrium properties. Let $\Delta$ be the overall $6\times 6$ stochastic amplitude matrix entering Equation (11),

$$\Delta =\left(\begin{array}{ccc}\frac{\mathsf{\alpha }}{m}& & \frac{\gamma }{m}\\ \frac{\mathsf{\epsilon }}{I}& & \frac{\beta }{I}\end{array}\right)$$

(12)

and define the $6\times 6$ matrix $\mathsf{\sigma }$, the entries of which are ${\sigma }_{i,j}$ as

$${\sigma }_{i,j}=\frac{1}{2}\sum _{h=1}^6{\Delta }_{i,h}{\Delta }_{j,h}$$

(13)

In matrix form,

$$\begin{array}{c}\hfill \mathsf{\sigma }=\frac{1}{2}\Delta {\Delta }^T=\left(\begin{array}{ccc}\frac{\mathsf{\alpha }{\mathsf{\alpha }}^T+\gamma {\gamma }^T}{2m^2}& & \frac{\alpha {\epsilon }^T+\gamma {\beta }^T}{2mI}\\ \frac{\epsilon {\alpha }^T+\beta {\gamma }^T}{2mI}& & \frac{\beta {\beta }^T+\epsilon {\epsilon }^T}{2I^2}\end{array}\right)=\left(\begin{array}{ccc}\frac{\mathbf{a}}{2m^2}& & \frac{\mathbf{c}}{2mI}\\ \frac{\mathbf{d}}{2mI}& & \frac{\mathbf{b}}{2I^2}\end{array}\right)\end{array}$$

(14)

which is, by definition, symmetric. Expressing Equation (11) componentwise

$$\begin{array}{ccc}\hfill d{v}_i& =& -\left(\frac{1}{m}\sum _{j=1}^3{\eta }_{i,j}{v}_j+\frac{1}{m}\sum _{j=1}^3{C}_{i,j}^{\left(1\right)}{\omega }_j\right)+\frac{1}{m}\sum _{j=1}^3{\alpha }_{i,j}d{w}_j+\frac{1}{m}\sum _{j=1}^3{\gamma }_{i,j}d{w}_j^{\omega }\hfill \\ \hfill d{\omega }_i& =& -\left(\frac{1}{I}\sum _{j=1}^3{C}_{i,j}^{\left(2\right)}{v}_j+\frac{1}{I}\sum _{j=1}^3{\eta }_{i,j}^{\omega }{\omega }_j\right)+\frac{1}{I}\sum _{j=1}^3{\epsilon }_{i,j}d{w}_j+\frac{1}{m}\sum _{j=1}^3{\beta }_{i,j}d{w}_j^{\omega }\hfill \end{array}$$

(15)

so that the associated Fokker-Planck equation for the probability density function $p(\mathbf{v},\mathsf{\omega },t)$ attains the form [45]

$$\begin{array}{ccc}\hfill \frac{\partial p}{\partial t}& =& \sum _{i=1}^3\frac{\partial }{\partial v_i}\left[\left(\frac{1}{m}\sum _{j=1}^3{\eta }_{i,j}{v}_j+\frac{1}{m}\sum _{j=1}^3{C}_{i,j}^{\left(1\right)}{\omega }_j\right)p\right]\hfill \\ \hfill & +& \sum _{i=1}^3\frac{\partial }{\partial {\omega }_i}\left[\left(\frac{1}{I}\sum _{j=1}^3{C}_{i,j}^{\left(2\right)}{v}_j+\frac{1}{I}\sum _{j=1}^3{\eta }_{i,j}^{\omega }{\omega }_j\right)p\right]\hfill \\ \hfill & +& \sum _{i,j=1}^3\frac{{\partial }^2}{\partial v_i\partial v_j}\left(\frac{a_{i,j}p}{2m^2}\right)+\sum _{i,j=1}^3\frac{{\partial }^2}{\partial v_i\partial {\omega }_j}\left(\frac{c_{i,j}p}{2mI}\right)\hfill \\ \hfill & +& \sum _{i,j=1}^3\frac{{\partial }^2}{\partial {\omega }_i\partial v_j}\left(\frac{d_{i,j}p}{2mI}\right)+\sum _{i,j=1}^3\frac{{\partial }^2}{\partial {\omega }_i\partial {\omega }_j}\left(\frac{b_{i,j}p}{2I^2}\right)\hfill \end{array}$$

(16)

The statistical equilibrium properties can be ascertained from the analysis of the lower-order (first and second) moments of $p(\mathbf{v},\mathsf{\omega },t)$. The first-order moments relax to zero in the long-term regime, as a consequence of the fact that the hydrodynamic interaction matrix $\mathbf{H}$ is positive definite. It is therefore sufficient to consider the second-order moments for the translational/angular velocities,

$$\begin{array}{ccc}\hfill M_{i,j}^{v,v}\left(t\right)& =& {\int }_{{\mathbb{R}}^3}d\mathsf{\omega }{\int }_{{\mathbb{R}}^3}{v}_i{v}_{j}p(\mathbf{v},\mathsf{\omega },t)d\mathbf{v}\hfill \\ \hfill M_{i,j}^{v,\omega }\left(t\right)& =& {\int }_{{\mathbb{R}}^3}{\omega }_{j}d\mathsf{\omega }{\int }_{{\mathbb{R}}^3}{v}_{i}p(\mathbf{v},\mathsf{\omega },t)d\mathbf{v}\hfill \\ \hfill M_{i,j}^{\omega ,\omega }\left(t\right)& =& {\int }_{{\mathbb{R}}^3}d\mathbf{v}{\int }_{{\mathbb{R}}^3}{\omega }_i{\omega }_{j}p(\mathbf{v},\mathsf{\omega },t)d\mathsf{\omega }\hfill \end{array}$$

(17)

To begin with, consider ${\mathbf{M}}^{v,v}$. From the Fokker-Planck Equation (16) one obtains

$$\begin{array}{ccc}\hfill \frac{d{M}_{h,k}^{v,v}}{dt}& =& -\frac{1}{m}\sum _{j=1}^3{\eta }_{h,j}{M}_{j,k}^{v,v}-\frac{1}{m}\sum _{j=1}^3{\eta }_{k,j}{M}_{j,h}^{v,v}-\frac{1}{m}\sum _{j=1}^3{C}_{h,j}^{\left(1\right)}{M}_{k,j}^{v,\omega }-\frac{1}{m}\sum _{j=1}^3{C}_{k,j}^{\left(1\right)}{M}_{h,j}^{v,\omega }\hfill \\ \hfill & +& \frac{a_{h,k}+a_{k,h}}{2m^2}\hfill \end{array}$$

(18)

In the long-term limit (equilibrium), it follows from Equation (18) that

$$\left[\mathsf{\eta }{\mathbf{M}}^{v,v}+{\mathbf{M}}^{v,v}\mathsf{\eta }\right]+\left[{\mathbf{C}}^{\left(1\right)}{\left({\mathbf{M}}^{v,\omega }\right)}^T+{\mathbf{M}}^{v,\omega }{\left({\mathbf{C}}^{\left(1\right)}\right)}^T\right]=\frac{1}{2m}\left(\mathbf{a}+{\mathbf{a}}^T\right)$$

(19)

Next, consider ${\mathbf{M}}^{\omega ,\omega }$, the entries of which satisfy the equations

$$\begin{array}{ccc}\hfill \frac{d{M}_{h,k}^{\omega ,\omega }}{dt}& =& -\frac{1}{I}\sum _{j=1}^3{C}_{h,j}^{\left(2\right)}{M}_{j,k}^{v,\omega }-\frac{1}{I}\sum _{j=1}^3{C}_{k,j}^{\left(2\right)}{M}_{j,h}^{v,\omega }-\frac{1}{I}\sum _{j=1}^3{\eta }_{h,j}^{\omega }{M}_{j,k}^{\omega ,\omega }-\frac{1}{I}\sum _{j=1}^3{\eta }_{k,j}^{\omega }{M}_{j,h}^{\omega ,\omega }\hfill \\ \hfill & +& \frac{b_{h,k}+b_{k,h}}{2I^2}\hfill \end{array}$$

(20)

so that the value attained at equilibrium is

$$\left[{\mathsf{\eta }}^{\omega }{\mathbf{M}}^{\omega ,\omega }+{\mathbf{M}}^{\omega ,\omega }{\mathsf{\eta }}^{\omega }\right]+\left[{\mathbf{C}}^{\left(2\right)}{\mathbf{M}}^{v,\omega }+{\left({\mathbf{M}}^{v,\omega }\right)}^T{\left({\mathbf{C}}^{\left(2\right)}\right)}^T\right]=\frac{1}{2I}\left(\mathbf{b}+{\mathbf{b}}^T\right)$$

(21)

Finally, consider the mixed second-order roto-translational moments

$$\begin{array}{ccc}\hfill \frac{d{M}_{h,k}^{v,\omega }}{dt}& =& -\frac{1}{m}\sum _{j=1}^3{\eta }_{h,j}{M}_{j,k}^{v,\omega }-\frac{1}{m}{C}_{h,j}^{\left(1\right)}{M}_{j,k}^{\omega ,\omega }-\frac{1}{I}\sum _{j=1}^3{\eta }_{k,j}^{\omega }{M}_{h,j}^{v,\omega }-\frac{1}{I}\sum _{j=1}^3{C}_{k,j}^{\left(2\right)}{M}_{j,h}^{v,v}\hfill \\ \hfill & +& \frac{c_{h,k}+d_{k,h}}{2mI}\hfill \end{array}$$

(22)

admitting the equilibrium condition

$$\frac{1}{m}\left[\mathsf{\eta }{\mathbf{M}}^{v,\omega }+{\mathbf{C}}^{\left(1\right)}{\mathbf{M}}^{\omega ,\omega }\right]+\frac{1}{I}\left[{\mathbf{M}}^{v,\omega }{\mathsf{\eta }}^{\omega }+{\mathbf{M}}^{v,v}{\left({\mathbf{C}}^{\left(2\right)}\right)}^T\right]\frac{\mathbf{c}+{\mathbf{d}}^T}{2mI}$$

(23)

Fluctuation-dissipation conditions, and the explicit expression for the stochastic amplitude matrices follow by enforcing the equilibrium properties [14]

$${\mathbf{M}}^{v,v}=\frac{k_{B}T}{m}\mathbf{I},{\mathbf{M}}^{\omega ,\omega }=\frac{k_{B}T}{I}\mathbf{I},{\mathbf{M}}^{v,\omega }=0$$

(24)

where I is the moment of inertia, and $\mathbf{I}$ the identity matrix. These conditions stem from the Maxwellian equilibrium velocity distribution, and from the energy equipartition theorem applied to a rigid particle, admitting 6 degrees of freedom. Making use of Equation (24) and of the equilibrium expression Equation (19) for ${\mathbf{M}}^{v,v}$, one obtains that the matrix $\mathbf{a}$ should be symmetric and

$$\mathbf{a}=2k_{B}T\mathsf{\eta }$$

(25)

A similar analysis for ${\mathbf{M}}^{\omega ,\omega }$ Equation (21) at equilibrium provides

$$\mathbf{b}=2k_{B}T{\mathsf{\eta }}^{\omega }$$

(26)

The equilibrium results deriving from the analysis of the mixed roto-translational moments yield

$$k_{B}T\left[{\mathbf{C}}^{\left(1\right)}+{\left({\mathbf{C}}^{\left(2\right)}\right)}^T\right]=\frac{\mathbf{c}+{\mathbf{d}}^T}{2}$$

(27)

a solution of which is ${\mathbf{d}}^T=\mathbf{c}$ and

$$\mathbf{c}=2k_{B}T{\mathbf{C}}^{\left(1\right)},\mathbf{d}=2k_{B}T{\mathbf{C}}^{\left(2\right)}$$

(28)

Once the “diffusional” matrices $\mathbf{a}$, $\mathbf{b}$, $\mathbf{c}$, $\mathbf{d}$ have been expressed in terms of the hydrodynamic resistance matrices, the next step is to derive the expression for the stochastic amplitude matrices $\mathsf{\alpha }$, $\mathsf{\beta }$, $\mathsf{\gamma }$, $\mathsf{\epsilon }$.

Because of Equation (28), the $6\times 6$ matrix $\Delta$ defined by Equation (12) is symmetric, and satisfies the algebraic matrix equation

$${\displaystyle \frac{{\Delta }^2}{2}}=k_{B}T\left(\begin{array}{ccc}\frac{\mathsf{\eta }}{m^2}& & \frac{{\mathbf{C}}^{\left(1\right)}}{mI}\\ \frac{{\mathbf{C}}^{\left(2\right)}}{mI}& & \frac{{\mathsf{\eta }}^{\omega }}{I^2}\end{array}\right)=\mathsf{\sigma }$$

(29)

Mutuating this property from the symmetry of hydrodynamic matrices, the matrix $\mathsf{\sigma }$ is symmetric and positive definite, and it is known from matrix theory [46] that there exists a unique, symmetric, and positive definite matrix $\Delta$ solution of Equation (29), formally

$$\Delta =\sqrt{2}{\mathsf{\sigma }}^{1/2}$$

(30)

In order to determine the explicit expression for the matrix $\Delta$, it is convenient to normalize its entries, expressing the force/torque and the velocity/angular velocity in the same physical dimensions. If ${\ell }_p$ is the characteristic particle length, ${\ell }_p=R_p$ for spherical particle of radius $R_p$, set $\widehat{\mathbf{T}}=\mathbf{T}/{\ell }_p$, $\hat{\mathsf{\omega }}={\ell }_p\mathsf{\omega }$. In this way, $\widehat{\mathbf{T}}$ has the dimension of a force, and $\hat{\mathsf{\omega }}$ the dimension of a velocity, so that Equations (7) and (8) become

$$\left(\begin{array}{c}\mathbf{F}\\ \widehat{\mathbf{T}}\end{array}\right)=-\widehat{\mathbf{H}}\left(\mathbf{x}\right)\left(\begin{array}{c}\mathbf{v}\\ \widehat{\mathsf{\omega }}\end{array}\right)$$

(31)

where $\widehat{\mathbf{H}}\left(\mathbf{x}\right)$ is the normalized overall resistance matrix possessing the block structure

$$\widehat{\mathbf{H}}\left(\mathbf{x}\right)=\left(\begin{array}{ccc}\widehat{\mathsf{\eta }}\left(\mathbf{x}\right)& & {\widehat{\mathbf{C}}}^{\left(1\right)}\left(\mathbf{x}\right)\\ {\widehat{\mathbf{C}}}^{\left(2\right)}\left(\mathbf{x}\right)& & {\widehat{\mathsf{\eta }}}^{\omega }\left(\mathbf{x}\right)\end{array}\right)=\left(\begin{array}{ccc}\mathsf{\eta }\left(\mathbf{x}\right)& & {\mathbf{C}}^{\left(1\right)}\left(\mathbf{x}\right)/{\ell }_p\\ {\mathbf{C}}^{\left(2\right)}\left(\mathbf{x}\right)/{\ell }_p& & {\mathsf{\eta }}^{\omega }\left(\mathbf{x}\right)/{\ell }_p^2\end{array}\right)$$

(32)

and the normalized moment of inertia is given by $\widehat{I}=I/{\ell }_p^2$. In this way, Equation (11) attains the normalized representation

$$\begin{array}{ccc}\hfill d\mathbf{v}\left(t\right)& =& -\frac{\eta }{m}\mathbf{v}\left(t\right)dt-\frac{{\widehat{\mathbf{C}}}^{\left(1\right)}}{m}\widehat{\mathsf{\omega }}\left(t\right)dt+\frac{\widehat{\alpha }}{m}d\mathbf{w}\left(t\right)+\frac{\widehat{\gamma }}{m}d{\mathbf{w}}^{\omega }\left(t\right)\hfill \\ \hfill d\widehat{\mathsf{\omega }}\left(t\right)& =& -\frac{{\widehat{\mathbf{C}}}^{\left(2\right)}}{\widehat{I}}\mathbf{v}\left(t\right)dt-\frac{{\widehat{\eta }}^{\omega }}{\widehat{I}}\widehat{\mathsf{\omega }}\left(t\right)dt+\frac{\widehat{\epsilon }}{\widehat{I}}d\mathbf{w}\left(t\right)+\frac{\widehat{\beta }}{\widehat{I}}d{\mathbf{w}}^{\omega }\left(t\right)\hfill \end{array}$$

(33)

with $\hat{\mathsf{\alpha }}=\mathsf{\alpha }$, $\hat{\mathsf{\gamma }}=\gamma /{\ell }_p$, $\hat{\mathsf{\epsilon }}=\mathsf{\epsilon }/{\ell }_p$, $\hat{\beta }=\beta /{\ell }_p^2$. The normalized matrix Equation (29) thus becomes

$${\displaystyle \frac{{\hat{\Delta }}^2}{2}}=k_{B}T\left(\begin{array}{ccc}\frac{\mathsf{\eta }}{{\mathbf{m}}^2}& & \frac{{\widehat{\mathbf{C}}}^{\left(1\right)}}{m\widehat{I}}\\ \frac{{\widehat{\mathbf{C}}}^{\left(2\right)}}{m\widehat{I}}& & \frac{{\widehat{\mathsf{\eta }}}^{\mathsf{\omega }}}{{\widehat{I}}^2}\end{array}\right)=\widehat{\mathsf{\sigma }}$$

(34)

and $\hat{\Delta }=\sqrt{2}{\widehat{\sigma }}^{1/2}$.

For symmetric and positive definite matrices $\hat{\mathsf{\sigma }}$, the estimate of their square root ${\hat{\mathsf{\sigma }}}^{1/2}$ reduces to an eigenvalue problem [46]. Let ${\lambda }_i>0$, $i=1,\cdots ,6$, be the eigenvalues of $\widehat{\mathsf{\sigma }}$, and ${\mathbf{V}}^{\left(i\right)}=(V_1^{\left(i\right)},\cdots ,V_6^{\left(i\right)})$ the corresponding unit eigenvectors. The solution of Equation (34) can be expressed as

$$\hat{\Delta }=\sqrt{2}\mathbf{V}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\lambda }_1^{1/2},\cdots ,{\lambda }_6^{1/2}){\mathbf{V}}^{-1}$$

(35)

where $\mathbf{V}$ is the eigenbasis transformation matrix, the column of which are orderly the eigevectors of $\widehat{\sigma }$, i.e.,

$$\mathbf{V}=\left(\begin{array}{ccc}{V}_1^{\left(1\right)}& \cdots & V_1^{\left(6\right)}\\ \cdots & \cdots & \cdots \\ V_6^{\left(1\right)}& \cdots & V_6^{\left(6\right)}\end{array}\right)$$

(36)

and $\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\lambda }_1^{1/2},\cdots ,{\lambda }_6^{1/2})$ is a diagonal matrix, the diagonal entries of which are the square roots of the eigenvalues ${\lambda }_i$, $i=1,\cdots ,6$.

An Example: Systems with Axialsymmetric Geometry

As an application of the previous analysis, consider a system with axialsymmetric geometry. Typical axisymmetric systems in microfluidics are spherical particles moving in a slit channel or near an infinitely extended planar wall. The symmetries of the problem reduces the 21 independent coefficients of the hydrodynamic resistance matrix to 5. By taking a Cartesian coordinate system $(x_1,x_2,x_3)$ with $x_3$ lying on the axis of symmetry as in Figure 1, the resistance matrix takes the form [33]

$$\mathbf{H}=\left(\begin{array}{cccccc}{\eta }_{1,1}& 0& 0& 0& C& 0\\ 0& {\eta }_{1,1}& 0& -C& 0& 0\\ 0& 0& {\eta }_{3,3}& 0& 0& 0\\ 0& -C& 0& {\eta }_{1,1}^{\omega }& 0& 0\\ C& 0& 0& 0& {\eta }_{1,1}^{\omega }& 0\\ 0& 0& 0& 0& 0& {\eta }_{3,3}^{\omega }\end{array}\right)$$

(37)

Therefore, the matrix $\hat{\mathsf{\sigma }}$, chosing ${\ell }_p=R_p$, becomes

$$\widehat{\mathsf{\sigma }}=\left(\begin{array}{cccccc}{\widehat{\sigma }}_1& 0& 0& 0& {\widehat{\sigma }}_c& 0\\ 0& {\widehat{\sigma }}_1& 0& -{\widehat{\sigma }}_c& 0& 0\\ 0& 0& {\widehat{\sigma }}_3& 0& 0& 0\\ 0& -{\widehat{\sigma }}_c& 0& {\widehat{\sigma }}_4& 0& 0\\ {\widehat{\sigma }}_c& 0& 0& 0& {\widehat{\sigma }}_4& 0\\ 0& 0& 0& 0& 0& {\widehat{\sigma }}_6\end{array}\right)$$

(38)

where

$$\begin{array}{ccc}\hfill {\widehat{\sigma }}_1& =& {\widehat{\sigma }}_2={\eta }_{1,1}{\displaystyle \frac{k_{B}T}{m^2}},{\widehat{\sigma }}_3={\eta }_{3,3}{\displaystyle \frac{k_{B}T}{m^2}},{\widehat{\sigma }}_4={\widehat{\sigma }}_5={\eta }_{1,1}^{\omega }{\displaystyle \frac{k_{B}T}{I^2}}{R}_p^2\hfill \\ \hfill {\widehat{\sigma }}_6& =& {\eta }_{3,3}^{\omega }{\displaystyle \frac{k_{B}T}{I^2}}{R}_p^2,{\widehat{\sigma }}_c=C{\displaystyle \frac{k_{B}T}{mI}}{R}_p\hfill \end{array}$$

(39)

The eigenvalues of $\hat{\mathsf{\sigma }}$ are

$$\begin{array}{c}\hfill {\lambda }_1={\lambda }_2={\displaystyle \frac{{\widehat{\sigma }}_1+{\widehat{\sigma }}_4-r}{2}},{\lambda }_3={\widehat{\sigma }}_3,{\lambda }_4={\lambda }_5={\displaystyle \frac{{\widehat{\sigma }}_1+{\widehat{\sigma }}_4+r}{2}},{\lambda }_6={\widehat{\sigma }}_6\end{array}$$

(40)

where $r=\sqrt{{({\widehat{\sigma }}_4-{\widehat{\sigma }}_1)}^2+4{\widehat{\sigma }}_c^2}$. The eigenvector matrix $\mathbf{V}$ entering Equation (36) takes in this case the expression

$$\mathbf{V}=\left(\begin{array}{cccccc}-\frac{-{\widehat{\sigma }}_1+{\widehat{\sigma }}_4+r}{2{\widehat{\sigma }}_c}& 0& 0& -\frac{-{\widehat{\sigma }}_1+{\widehat{\sigma }}_4-r}{2{\widehat{\sigma }}_c}& 0& 0\\ 0& -\frac{{\widehat{\sigma }}_1-{\widehat{\sigma }}_4-r}{2{\widehat{\sigma }}_c}& 0& 0& -\frac{{\widehat{\sigma }}_1-{\widehat{\sigma }}_4+r}{2{\widehat{\sigma }}_c}& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 1& 0& 0& 1& 0\\ 1& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right)$$

(41)

By applying Equations (36) and (38), we obtain for $\hat{\Delta }$

$$\hat{\Delta }=\left(\begin{array}{cccccc}{\hat{\Delta }}_1& 0& 0& 0& {\hat{\Delta }}_c& 0\\ 0& {\hat{\Delta }}_1& 0& -{\hat{\Delta }}_c& 0& 0\\ 0& 0& {\hat{\Delta }}_3& 0& 0& 0\\ 0& -{\hat{\Delta }}_c& 0& {\hat{\Delta }}_4& 0& 0\\ {\hat{\Delta }}_c& 0& 0& 0& {\hat{\Delta }}_4& 0\\ 0& 0& 0& 0& 0& {\hat{\Delta }}_6\end{array}\right)$$

(42)

where

$$\begin{array}{ccc}\hfill {\hat{\Delta }}_1& =& {\hat{\Delta }}_2={\displaystyle \frac{({\widehat{\sigma }}_1-{\lambda }_1)\sqrt{{\lambda }_4}-({\widehat{\sigma }}_1-{\lambda }_4)\sqrt{{\lambda }_1}}{r}},{\widehat{\Delta }}_3=\sqrt{2{\widehat{\sigma }}_3}\hfill \\ \hfill {\widehat{\Delta }}_4& =& {\widehat{\Delta }}_5={\displaystyle \frac{({\widehat{\sigma }}_4-{\lambda }_4)\sqrt{{\lambda }_1}-({\widehat{\sigma }}_4-{\lambda }_1)\sqrt{{\lambda }_4}}{r}},{\widehat{\Delta }}_6=\sqrt{2{\widehat{\sigma }}_6},{\widehat{\Delta }}_c={\displaystyle \frac{{\widehat{\sigma }}_c(\sqrt{{\lambda }_4}-\sqrt{{\lambda }_1})}{r}}\hfill \end{array}$$

(43)

As $h\to \infty$, i.e., far away from the wall, the coupling term ${\hat{\mathsf{\sigma }}}_c\to 0$, and the quantity r take the limit form

$$r=\sqrt{{({\widehat{\sigma }}_4-{\widehat{\sigma }}_1)}^2}={\displaystyle \frac{\pi \mu R_p{k}_{B}T}{m^2}}\sqrt{{\left({\displaystyle \frac{{\overline{\eta }}_{1,1}^{\omega }}{{\overline{I}}^2}}-{\overline{\eta }}_{1,1}\right)}^2}$$

(44)

where

$${\overline{\eta }}_{1,1}^{\omega }={\displaystyle \frac{{\eta }_{11}^{\omega }}{\pi \mu R_p^3}}=8,{\overline{\eta }}_{1,1}={\displaystyle \frac{{\eta }_{11}}{\pi \mu R_p}}=6,\overline{I}={\displaystyle \frac{I}{m{R}_p^2}}={\displaystyle \frac{2}{5}}$$

(45)

Substituting these values into Equation (43), we obtain the expected results in the free space

$${\widehat{\Delta }}_i=\sqrt{2{\widehat{\sigma }}_i},i=1,\cdots ,6$$

(46)

and therefore

$$\mathsf{\alpha }=\sqrt{2k_{B}T{\eta }_{\infty }}\mathbf{I},\mathsf{\beta }=\sqrt{2k_{B}T{\eta }_{\infty }^{\omega }}\mathbf{I},\mathsf{\gamma }=\mathsf{\epsilon }=0$$

(47)

where $\mathbf{I}$ is the $3\times 3$ identity matrix. In order to recover the Stokesian limit values in the free space ${\eta }_{\infty }$, ${\eta }_{\infty }^{\omega }$, the argument of the square root in Equation (43) should be positive, i.e.,

$${\displaystyle \frac{{\overline{\eta }}_{1,1}^{\omega }}{{\overline{\eta }}_{1,1}}}>{\overline{I}}^2$$

(48)

Inequality (48) is fulfilled in the free space as follows from Equation (45).

3. Adiabatic Elimination of the Velocity Variables

In this section we consider the formulation of the overdamped approximation for micrometric rigid particles in confined geometries. The overdamped approximation consists in expressing the equations of motion in the form of a kinematic equation for the mechanical degrees of freedom associated with translational and rotational motions. This is made possible due to the fact that velocity variables are customarily characterized by a faster relaxation dynamics than position and orientational variables. This is certainly true for micrometric particles if one considers their transport properties at time scales much larger than the characteristic dissipation timescale $t_{\mathrm{diss}}=m/{\eta }_{\infty }$, the order of magnitude of which falls between ${10}^{-6}$–${10}^{-7}$ s for micrometric particles in water at room temperature. For this reason, the overdamped approximation is often referred to as the adiabatic elimination of the fast velocity variables, and this follows by imposing the condition

$$md\mathbf{v}\simeq 0,Id\mathsf{\omega }\simeq 0$$

(49)

and extracting out of Equation (49) the expression for $\mathbf{v}$ and $\mathsf{\omega }$, entering the kinematic equation

$$d\mathbf{x}=\mathbf{v}dt,d\mathsf{\varphi }=\omega dt$$

(50)

where $\mathsf{\varphi }=({\varphi }_1,{\varphi }_2,{\varphi }_3)$ is the vector-valued angular variable accounting for the particle orientation.

The overdamped approximation in confined geometries presents intrinsic peculiarities, just because the hydrodynamic resistance matrix depends on the position $\mathbf{x}$, and in general of the orientation $\mathsf{\varphi }$, and this raises delicate issues when the thermal fluctuations are expressed as linear superposition of increments of Wiener processes, owing to their highly singular local structure [45]. This problem in the physical literature is usually referred to as the Ito-Stratonovich dilemma [47,48]. The most convenient and simple approach to perform the adiabatic elimination of the fast velocity variable is due to Sancho et al. [49]. The starting point in this derivation is that the configurational coordinates of a particle driven by Wiener fluctuations still represent a locally smooth, and almost everywhere differentiable continuous stochastic process with probability 1, and this property determines the way Equation (49) is interpreted. Without loss of generality, let us suppose that the particle is subjected to an external potential $U\left(\mathbf{x}\right)$, and that no additional flow contribution are present. The latter can be added at the end of the adiabatic elimination process.

In order to perform this analysis in the simplest formal way, it is convenient to group together configurational and velocity variables, thus introducing the 6-dimensional configurational and velocity variables, $\mathbf{z}$ and $\mathbf{V}$, respectively, and the overall stochastic forcing $d\mathbf{W}\left(t\right)$

$$\mathbf{z}=\left(\begin{array}{c}\mathbf{x}\\ \mathsf{\varphi }\end{array}\right),\mathbf{V}=\left(\begin{array}{c}\mathbf{v}\\ \mathsf{\omega }\end{array}\right),d\mathbf{W}\left(t\right)=\left(\begin{array}{c}d\mathbf{w}\left(t\right)\\ d{\mathbf{w}}^{\omega }\left(t\right)\end{array}\right)$$

(51)

so that the equations of motion can be compactly expressed as

$$d\mathbf{z}=\mathbf{V}dt,{\mathbf{m}}_{mI}d\mathbf{z}=-\mathbf{H}\left(\mathbf{z}\right)\mathbf{V}dt-{\nabla }_{z}U\left(\mathbf{z}\right)dt+\nu \left(\mathbf{z}\right)d\mathbf{W}\left(t\right)$$

(52)

where ${\mathbf{m}}_{mI}$ is the mass/moment-of inertial tensor, ${\mathbf{m}}_{mI}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(m,m,m,I,I,I)$ for a spherical particle, $\mathbf{H}\left(\mathbf{z}\right)$ the hydrodynamic resistance matrix, ${\nabla }_z$ the nabla operator with respect to the $\mathbf{z}$-variable and $\mathsf{\nu }$ is the matrix of thermal fluctuation intensity satisfying, as discussed in the previous Section, the condition at thermal equilibrium,

$$\mathsf{\nu }\left(\mathbf{z}\right){\mathsf{\nu }}^T\left(\mathbf{z}\right)=2k_{B}T\mathbf{H}\left(\mathbf{z}\right)$$

(53)

The overdamped approximation corresponds to the limit for mass and momentum of inertia tending to zero. In performing this limit, upon correct physical grounds, it should be ensured that the configurational variable $\mathbf{z}\left(t\right)$ is a smooth stochastic process. To this end, the Wong-Zakai theorem can be enforced [50,51], implying that the stochastic differential equation describing particle dynamics should be interpreted in a Stratonovich way [52], Technically, this means that, in performing the limit process, the quantity $\mathbf{H}\left(\mathbf{z}\right)\mathbf{V}dt$ should be interpreted as

$$\mathbf{H}\left(\mathbf{z}\right)\mathbf{V}dt=\mathbf{H}\left(\mathbf{z}\right)\circ d\mathbf{z}=\mathbf{H}(\mathbf{z}+d\mathbf{z}/2)d\mathbf{z}$$

(54)

where “∘” indicates the Stratonovich rule in the definition of stochastic integrals and differentials. It is also convenient to recall a known result, deriving from the Ito lemma, namely that [45]

$$d{W}_i\left(t\right)d{W}_j\left(t\right)={\delta }_{i,j}dt+o\left(dt\right)$$

(55)

where ${\delta }_{i,j}$ are the Kronecker’s symbols (entries of the identity matrix), and $o\left(dt\right)$ is a quantity going to zero for $dt\to 0$ faster than $dt$.

The final goal of this analysis is to derive a kinematic equation for the configurational particle degrees of freedom (Langevin equation) expressed in the Ito way and, out of it, the transport equation for particle density, corresponding to the Fokker-Planck equation for the statistical characterization of the so-obtained Langevin equation.

Making use of the Wong-Zakai result, expressed by Equation (54), the overdamped approximation of Equation (52) is thus given componentwise by

$$0=\sum _{j=1}^6{H}_{i,j}\left(\mathbf{z}+\frac{\mathbf{dz}}{2}\right)d{z}_j-\frac{\partial U}{\partial z_i}dt+\sum _{j=1}^6{\nu }_{i,j}\left(\mathbf{z}\right)d{W}_j\left(t\right)$$

(56)

Expanding the first term at the r.h.s. of Equation (56) in Taylor series to the leading order,

$$\sum _{j=1}^6{H}_{i,j}\left(\mathbf{z}+\frac{\mathbf{dz}}{2}\right)d{z}_j=\sum _{j=1}^6{H}_{i,j}\left(\mathbf{z}\right)d{z}_j+\frac{1}{2}\sum _{j,k=1}^6\frac{\partial H_{i,j}\left(\mathbf{z}\right)}{\partial z_k}d{z}_{j}d{z}_k+o\left(dt\right)$$

(57)

The quantity $d{z}_{j}d{z}_k$ can be evaluated from Equation (56), by considering for the first term at the r.h.s. of Equation (56) its Ito interpretation, namely ${\sum }_{j=1}^6{H}_{i,j}\left(\mathbf{z}\right)d{z}_j$ that provides the leading order contribution as the remainder in this approximation is order of $o\left(dt\right)$. Thus, enforcing also Equation (55), one obtains

$$\begin{array}{ccc}\hfill d{z}_{j}d{z}_k& =& \left(-\sum _{p=1}^6{H}_{j,p}^{-1}\frac{\partial U}{\partial z_p}dt+\sum _{p,m=1}^6{H}_{j,p}^{-1}{\nu }_{p,m}d{W}_m\left(t\right)\right)\hfill \\ & \times & \left(-\sum _{q=1}^6{H}_{k,q}^{-1}\frac{\partial U}{\partial z_q}dt+\sum _{q,n}^6{H}_{k,q}^{-1}{\nu }_{q,n}d{W}_n\left(t\right)\right)\hfill \\ & =& \sum _{p,m,q,n=1}^6{H}_{j,p}^{-1}{H}_{k,q}^{-1}{\nu }_{p,m}{\nu }_{q,n}d{W}_m\left(t\right)d{W}_n\left(t\right)\hfill \\ & =& \sum _{p,m,q,n=1}^6{H}_{j,p}^{-1}{H}_{k,q}^{-1}{\nu }_{p,m}{\nu }_{q,n}{\delta }_{m,n}dt=\sum _{p,q,m=1}^6{H}_{j,p}^{-1}{H}_{k,q}^{-1}{\nu }_{p,m}{\nu }_{q,m}dt\hfill \\ & =& \sum _{p,q=1}^6{H}_{j,p}^{-1}{H}_{k,q}^{-1}{k}_{B}T{H}_{q,p}dt=2k_{B}T{H}_{j,k}^{-1}dt+o\left(dt\right)\hfill \end{array}$$

(58)

where $H_{j,p}^{-1}={\left({\mathbf{H}}^{-1}\right)}_{j,p}$, and in deriving the last relation we have made use of Equation (53). Substituting Equations (57) and (58) into Equation (56) it follows that

$$\sum _{j=1}^6{H}_{i,j}d{z}_j=-k_{B}T\sum _{j,k=1}^6\frac{\partial H_{i,j}}{\partial z_k}{H}_{k,j}^{-1}dt-\frac{\partial U}{\partial z_i}dt+\sum _{j=1}^6{\nu }_{i,j}d{W}_j\left(t\right)$$

(59)

where the $o\left(dt\right)$-terms have been neglected. Correspondingly, the Langevin equation in the configuration $(\mathbf{x},\mathsf{\varphi })$-space becomes

$$d{z}_i=-\sum _{j=1}^6{H}_{i,j}^{-1}\frac{\partial U}{\partial z_j}dt-f_{i}dt+\sum _{j,h=1}^6{H}_{i,j}^{-1}{\nu }_{j,h}d{W}_h\left(t\right)$$

(60)

where

$$f_i=k_{B}T\sum _{j,h,k=1}^6{H}_{i,j}^{-1}\frac{\partial H_{j,h}}{\partial z_k}{H}_{k,h}^{-1}$$

(61)

The Fokker-Planck equation for the probability density $p(\mathbf{z},t)$ associated with Equation (61) is thus given by

$$\frac{\partial p}{\partial t}={\nabla }_z·\left({\mathbf{H}}^{-1}{\nabla }_{z}Up\right)+{\nabla }_z·\left(\mathbf{f}p\right)+\sum _{i,j=1}^6\frac{{\partial }^2}{\partial z_i\partial z_j}\left(D_{i,j}p\right)$$

(62)

where $\mathbf{f}={\left(f_i\right)}_{i=1}^6$ and the generalized diffusivity tensor $D_{i,j}$ takes the form

$$\begin{array}{ccc}\hfill D_{i,j}& =& \frac{1}{2}\sum _{p,k,q=1}^6{H}_{i,p}^{-1}{\nu }_{p,k}{H}_{j,q}^{-1}{\nu }_{q,k}=\frac{1}{2}\sum _{p,q=1}^6{H}_{i,p}^{-1}{H}_{j,q}^{-1}2k_{B}T{H}_{p,q}\hfill \\ & =& k_{B}T{H}_{i,j}^{-1}\hfill \end{array}$$

(63)

i.e., the generalized diffusivity tensor $\mathbf{D}\left(\mathbf{z}\right)={\left(D_{i,j}\left(\mathbf{z}\right)\right)}_{i,j=1}^6$ is related to the resistance matrix $\mathbf{H}\left(\mathbf{z}\right)$ by the relation

$$\mathbf{D}\left(\mathbf{z}\right)\mathbf{H}\left(\mathbf{z}\right)=k_{B}T$$

(64)

generalizing the fluctuation-dissipation relation Equation (4). Next, consider the contribution of the vector field $\mathbf{f}$ entering the Fokker-Planck Equation (62). From the identity ${\mathbf{H}}^{-1}\mathbf{H}=\mathbf{I}$, it follows componentwise, for any $k=1,\cdots ,6$,

$$\sum _{j=1}^6\frac{\partial H_{i,j}^{-1}}{\partial z_k}{H}_{j,h}+\sum _{j=1}^6{H}_{i,j}^{-1}\frac{\partial H_{j,h}}{\partial z_k}=0$$

(65)

and thus,

$$\frac{\partial H_{j,h}^{-1}}{\partial z_k}=-\sum _{m,j=1}^6{H}_{i,m}^{-1}\frac{\partial H_{m,j}}{\partial z_k}{H}_{j,h}^{-1}$$

(66)

The latter expression implies that the entries $f_i$ of $\mathbf{f}$ defined by Equation (61) reduce to

$$f_i=-k_{B}T\sum _{h=1}^6\frac{\partial H_{i,h}^{-1}}{\partial z_h}$$

(67)

Substituting Equations (63) and (67) into Equation (62), the Fokker-Planck equation for $p(\mathbf{z},t)$ attains the simpler form

$$\frac{\partial p}{\partial t}=\sum _{i,j=1}^6\frac{\partial }{\partial z_i}\left(H_{i,j}^{-1}\frac{\partial U}{\partial z_j}p\right)-k_{B}T\sum _{i=1}^6\frac{\partial }{\partial z_i}\left(\sum _{j=1}^6\frac{\partial H_{i,j}^{-1}}{\partial z_j}p\right)+k_{B}T\sum _{i,j=1}^6\frac{\partial }{\partial z_i\partial z_j}\left(H_{i,j}^{-1}p\right)$$

(68)

that can be rewritten in a more compact way as

$$\frac{\partial p}{\partial t}=\sum _{i,j=1}^6\frac{\partial }{\partial z_i}\left(H_{i,j}^{-1}\frac{\partial U}{\partial z_j}p\right)+\sum _{i,j=1}^6\frac{\partial }{\partial z_i}\left(D_{i,j}\frac{\partial p}{\partial z_j}\right)$$

(69)

The latter corresponds to an advection-diffusion equation in the configuration space in the presence of the effective velocity ${\mathbf{v}}_{\mathrm{eff}}={\mathbf{H}}^{-1}{\nabla }_{z}U$, stemming from the potential $U\left(\mathbf{z}\right)$ and of the tensor diffusivity $\mathbf{D}$. Conversely, Equation (68) represents the classical formulation of the Fokker-Planck equation associated with a Langevin dynamics interpreted in the Ito way, attaining the following expression

$$d{z}_i=-\sum _{j=1}^6{H}_{i,j}^{-1}\frac{\partial U}{\partial z_j}dt+\sum _{j=1}^6\frac{\partial D_{i,h}}{\partial z_h}dt+\sqrt{2}\sum _{j=1}^6{\left({\mathbf{D}}^{1/2}\right)}_{i,j}d{W}_j\left(t\right)$$

(70)

where ${\left({\mathbf{D}}^{1/2}\right)}_{i,j}$ are the entries of the square root matrix ${\mathbf{D}}^{1/2}$ of the diffusivity tensor $\mathbf{D}$, ${\mathbf{D}}^{1/2}{\mathbf{D}}^{1/2}=\mathbf{D}$. It is also clear from the Ito representation of the Langevin Equation (70) the occurrence of an additional convective contribution depending on the divergence of the diffusivity tensor. This term admits a physical meaning as, even for $U\left(\mathbf{z}\right)=0$, it provides a biasing average velocity $V_i^{\left(\mathrm{bias}\right)}\left(\mathbf{z}\right)$,

$$V_i^{\left(\mathrm{bias}\right)}\left({\mathbf{z}}^{*}\right)={\left.\frac{d\langle z_i\rangle }{dt}\right|}_{\mathbf{z}={\mathbf{z}}^{*}}$$

(71)

where $d\langle z_i\rangle {/dt|}_{\mathbf{z}={\mathbf{z}}^{*}}$ is the average value of the particle velocity evaluated when the particle configuration is at $\mathbf{z}={\mathbf{z}}^{*}$ [53,54].

In the case of a spherical particle, the hydrodynamic matrices depends solely on $\mathbf{x}$ and not on $\mathsf{\varphi }$. Consequently it is easy to obtain the evolution equation for the marginal distribution $p_x(\mathbf{x},t)={\int }_{S_3}p(\mathbf{z},t)d\mathsf{\varphi }$, where $S_3={[0,2\pi )}^3$. Assume also that the spherical particles are immersed in a flow, and that the force exterted by the flow onto a generic particle located at $\mathbf{x}$ is ${\mathbf{F}}_{\mathrm{flow}}\left(\mathbf{x}\right)=\mathsf{\eta }\left(\mathbf{x}\right){\mathbf{u}}^{\left(p\right)}\left(\mathbf{x}\right)$. In this case, the spatial particle density function $p_x(\mathbf{x},t)$ satisfies the balance equation

$$\frac{\partial p_x}{\partial t}=\sum _{i=1}^3\frac{\partial }{\partial x_i}\left[\left(u_i^{\left(p\right)}-\sum _{j=1}^3{\eta }_{i,j}^{-1}\frac{\partial U}{\partial x_j}\right)p_x\right]+\sum _{i,j=1}^3\frac{\partial }{\partial x_i}\left(D_{i,j}^x\frac{\partial p_x}{\partial x_j}\right)$$

(72)

where $D_{i,j}^x\left(\mathbf{x}\right)$, is the $3\times 3$ diffusivity tensor, ${\sum }_{h=1}^3{D}_{i,h}^x\left(\mathbf{x}\right){\eta }_{h,j}\left(\mathbf{x}\right)=k_{B}T{\delta }_{i,j}$.

3.1. An Application

As a simple application of the overdamped theory, consider the vertical motion of a spherical particle in the upper half-plane delimited by a planar wall at $x_3=0$ in isothermal conditions at temperature T. Indicate with $x=h$ the distance of the particle from the wall, and assume that the particle is subjected to an external potential $U\left(x\right)$, as in [28,29], where $U\left(x\right)$ stems from gravity and from a local double-layer repulsive potential near the wall. In this case the problem is spatially one-dimensional, since $\eta \left(x\right)={\eta }_{3,3}\left(x\right)$, and $D\left(x\right)=k_{B}T/\eta \left(x\right)$ depends solely from the distance x from the wall. Setting $U^{\prime }\left(x\right)=dU\left(x\right)/dx$, and similarly for $D^{\prime }\left(x\right)$, the Fokker-Planck equation for the density function $p_x(x,t)$ reads

$$\frac{\partial p_x(x,t)}{\partial t}=\frac{\partial }{\partial x}\left(\frac{U^{\prime }\left(x\right)}{\eta \left(x\right)}{p}_x(x,t)\right)+\frac{\partial }{\partial x}\left(D\left(x\right)\frac{\partial p_x(x,t)}{\partial x}\right)$$

(73)

and this equation corresponds to the Langevin-Ito equation

$$dx\left(t\right)=-{\eta }^{-1}\left(x\left(t\right)\right)U^{\prime }\left(x\left(t\right)\right)dt+D^{\prime }\left(x\left(t\right)\right)dt+\sqrt{2D\left(x\right(t\left)\right)}dw\left(t\right)$$

(74)

where $dw\left(t\right)$ is the increment of a one-dimensional Wiener process. Assuming that the potential is attractive towards $x=0$ at large distances, so that ${\int }_0^{\infty }{e}^{-U\left(x\right)/k_{B}T}dx<\infty$, it follows from Equation (73) that in the limit for $t\to \infty$, the density $p_x(x,t)$ converges towards a stationary density $p_x^{*}\left(x\right)$, solution of the equation

$${\eta }^{-1}\left(x\right)U^{\prime }\left(x\right)p_x^{*}\left(x\right)+D\left(x\right)\frac{d{p}_x^{*}\left(x\right)}{dx}=0$$

(75)

corresponding, as expected, to the Boltzmann distribution

$$p_x^{*}\left(x\right)=A{e}^{-U\left(x\right)/k_{B}T}$$

(76)

This is a classical result, starting from which, micrometric particles may be used as Brownian probes, upon recording their statistical properties i.e. their stationary density function $p_x^{*}\left(x\right)$, in order to investigate and measure surface properties of materials [28,29].

It is instructive to analyze in greater detail the mathematical properties of Equation (73). This is a one-dimensional parabolic equation for $p_x(x,t)$, and as the wall at $x=0$ is impermeable to particle transport, ${\int }_0^{\infty }{p}_x(x,t)dx=1$ for any $t\ge 0$. To solve Equation (73), it should be equipped with initial and boundary conditions. As regards the initial condition, $p_x(x,0)=p_{x,0}\left(x\right)$, with ${\int }_0^{\infty }{p}_{x,0}\left(x\right)dx=1$. The boundary condition at infinity, $x\to \infty$ is the classical regularity condition, namely

$$\underset{x\to \infty }{lim}{x}^n{p}_x(x,t)=\underset{x\to \infty }{lim}{x}^n\frac{\partial p_x(x,t)}{\partial x}=0,\mathrm{f}\mathrm{o}\mathrm{r} t>0,n=0,1,\cdots$$

(77)

meaning that $p_x(x,t)$ and $\partial p_x(x,t)/\partial x$ should decay faster than any power $x^n$ of x for $x\to \infty$. At $x=0$, the zero-flux boundary condition applies, implying in the present case

$${\left.-\frac{U^{\prime }\left(x\right)}{\eta \left(x\right)}{p}_x(x,t)-D\left(x\right)\frac{\partial p_x(x,t)}{\partial x}\right|}_{x=0}=0$$

(78)

Two cases should be discussed. If, (i) ${lim}_{x\to 0}\left|U^{\prime }\left(x\right)\right|<\infty$, since $D\left(0\right)=k_{B}T/\eta \left(0\right)=0$, the wall boundary condition at $x=0$ blows up, reducing Equation (78) to a trivial identity that does not provide any condition on the local behavior of $p_x(x,t)$ near $x=0$. Conversely, if (ii) ${lim}_{x\to 0}\left|U^{\prime }\left(x\right)\right|=\infty$, and moreover ${lim}_{x\to 0}{U}^{\prime }\left(x\right)/\eta \left(x\right)=C\ne 0$, where the constant C may even diverge to ∞, the boundary condition (78) reduces to a homogeneous Dirichlet condition,

$$p_x(0,t)=0$$

(79)

As in principle, the condition on the potential U leading to Equation (79) could not be verified in a physical systems, as for the case analyzed in [28,29], it follows from the above analysis that the simple transport model Equation (73) displays a singular behavior as regards the wall boundary conditions. This singular phenomenon is a peculiar feature of the transport equations involving a hydrodynamic Stokesian description of the fluid-particle interactions in confined geometries, deriving from the singularity of some entries of the resistance matrix near a solid wall.

3.2. Thermophoresis from the Overdamped Approximation

An interesting byproduct of the overdamped analysis discussed above is the derivation of thermophoretic effects from the stochastic equations of motion. For simplicity, let us consider the case of a spherical particle in its translational motion, neglecting rotational effects and in the absence of any external or fluid forcing. It has been shown in [55] (see also [56]) that even in the presence of a non-equilibrium steady temperature profile $T\left(\mathbf{x}\right)$ the fluctuation-dissipation relation can be applied, so that the equations of motion in the present case read

$$\begin{array}{ccc}\hfill d\mathbf{x}& =& \mathbf{v}dt\hfill \\ \hfill md\mathbf{v}& =& -\mathsf{\eta }\left(\mathbf{x}\right)\mathbf{v}dt+\sqrt{2k_{B}T\left(\mathbf{x}\right)}{\mathsf{\eta }}^{1/2}\left(\mathbf{x}\right)d\mathbf{w}\left(t\right)\hfill \end{array}$$

(80)

Within the overdamped approximation, enforcing the same Stratonovich-approach developed in Section 2 to the term $\mathsf{\eta }\left(\mathbf{x}\right)\mathbf{v}dt=\mathsf{\eta }\left(\mathbf{x}\right)\circ d\mathbf{x}$, one obtains

$$-\sum _{j=1}^3{\eta }_{i,j}d{x}_j-\frac{1}{2}\sum _{j,k=1}^3\frac{\partial {\eta }_{i,j}}{\partial x_k}d{x}_{j}d{x}_k+\sqrt{2k_{B}T\left(\mathbf{x}\right)}{\left({\mathsf{\eta }}^{1/2}\right)}_{i,j}d{w}_j\left(t\right)=0$$

(81)

Since,

$$\begin{array}{ccc}\hfill d{x}_{j}d{x}_k& =& 2k_{B}T\left(\mathbf{x}\right)\sum _{p,q=1}^3{\left({\mathsf{\eta }}^{1/2}\right)}_{j,p}{\left({\mathsf{\eta }}^{1/2}\right)}_{k,q}d{w}_p\left(t\right)d{w}_q\left(t\right)+o\left(dt\right)\hfill \\ \hfill & =& 2k_{B}T\left(\mathbf{x}\right)\sum _{p,q=1}^3{\left({\mathsf{\eta }}^{1/2}\right)}_{j,p}{\left({\mathsf{\eta }}^{1/2}\right)}_{k,p}dt+o\left(dt\right)=2k_{B}T\left(\mathbf{x}\right){\eta }_{j,k}^{-1}dt+o\left(dt\right)\hfill \end{array}$$

(82)

This leads to the following Langevin equation

$$d{x}_i=-k_{B}T\left(\mathbf{x}\right)\sum _{j,h,k=1}^3{\eta }_{i,j}^{-1}\frac{\partial {\eta }_{j,h}}{\partial x_k}{\eta }_{h,k}^{-1}dt+\sqrt{2k_{B}T\left(\mathbf{x}\right)}{\left({\mathsf{\eta }}^{1/2}\right)}_{i,j}d{w}_j\left(t\right)$$

(83)

From the fluctuation-dissipation relation extended to non-equilibrium steady thermal conditions, ${\sum }_{h=1}^3{D}_{i,h}{\eta }_{h,j}=k_{B}T\left(\mathbf{x}\right){\delta }_{i,j}$, it follows, for any $k=1,2,3$, that

$$\begin{array}{ccc}\hfill \sum _{h=1}^3\frac{\partial D_{i,h}}{\partial x_k}{\eta }_{h,j}+\sum _{h=1}^3{D}_{i,h}\frac{\partial {\eta }_{h,j}}{\partial x_k}& =& k_B\frac{\partial T}{\partial x_k}{\delta }_{i,j}\hfill \\ \hfill \sum _{h=1}^3\frac{\partial D_{i,h}}{\partial x_k}{\eta }_{h,j}+k_{B}T\left(\mathbf{x}\right)\sum _{h=1}^3{\eta }_{i,h}^{-1}\frac{\partial {\eta }_{h,j}}{\partial x_k}& =& k_B\frac{\partial T}{\partial x_k}{\delta }_{i,j}\hfill \end{array}$$

(84)

which implies

$$-k_{B}T\left(\mathbf{x}\right)\sum _{h,j=1}^3{\eta }_{i,h}^{-1}\frac{\partial {\eta }_{h,j}}{\partial x_k}{\eta }_{j,m}=\frac{\partial D_{i,m}}{\partial x_k}-k_B\frac{\partial T}{\partial x_k}{\eta }_{i,m}^{-1}$$

(85)

From the latter expression it follows that

$$-k_{B}T\left(\mathbf{x}\right)\sum _{h,j,k=1}^3{\eta }_{i,h}^{-1}\frac{\partial {\eta }_{h,j}}{\partial x_k}{\eta }_{j,k}=\sum _{k=1}^3\frac{\partial D_{i,k}}{\partial x_k}-\frac{1}{T\left(\mathbf{x}\right)}\frac{\partial T}{\partial x_k}{D}_{i,k}$$

(86)

Therefore, Equation (82) becomes

$$d{x}_i=\left(\sum _{k=1}^3\frac{\partial D_{i,k}}{\partial x_k}\right)dt-\sum _{k=1}^3\frac{D_{i,k}}{T\left(\mathbf{x}\right)}\frac{\partial T}{\partial x_k}dt+\sum _{k=1}^3\sqrt{2}{\left({\mathbf{D}}^{1/2}\right)}_{i,k}d{w}_k\left(t\right)$$

(87)

and the corresponding Fokker-Planck equations reads

$$\frac{\partial p_x}{\partial t}=\sum _{i,j=1}^3\frac{\partial }{\partial x_i}\left(\frac{D_{i,j}}{T}\frac{\partial T}{\partial x_j}{p}_x\right)+\sum _{i,j=1}^3\frac{\partial }{\partial x_i}\left(D_{i,j}\frac{\partial p_x}{\partial x_j}\right)$$

(88)

providing the occurrence of an additional convective contribution to the flux, depending on the temperature gradient, and equal to

$${\mathbf{J}}_{\mathrm{thermo}}=-\frac{\mathbf{D}}{T}\nabla T{p}_x$$

(89)

providing a thermophoretic velocity ${\mathbf{v}}_{\mathrm{thermo}}$ equal to $-\mathbf{D}\nabla T/T$ [57]. This result shows that thermophoretic effects naturally follow from the accurate description of stochastic particle motion in thermal non-equilibrium conditions.

4. Wall Singularities: Superficial Phenomena

In the previous Section we have outlined the physical problems associated with the setting of boundary conditions for transport equations in confined geometries due to the singularities of the entries of the resistance matrix. The asymptotic trends of all the entries of $\mathbf{H}\left(\mathbf{x}\right)$ of two smooth no-slip surfaces almost in contact, such as the case of a particle very close to a solid surface, has been studied by Cox [58] for any value of the walls’ curvatures. From the Cox’s analysis, it results that the drag force on a particle moving towards the wall is inversely proportional to gap between the surfaces (particle and wall surfaces). To simplify the analysis, let us focus on a spherical particle of radius $R_p$. By considering a Cartesian reference system $\mathbf{x}=(x_1,x_2,x_3)$ with the origin on the wall and with $x_3$ colinear to the axis passing between the two contact points of the two surfaces, we have that ${\eta }_{3,3}\left(\mathbf{x}\right)\sim O(R_p/h)$ as $h/R_p\to 0$. For example, the drag force on a sphere moving towards a planar wall can be expressed to the leading order by the Taylor approximation [59]

$${\eta }_{3,3}\left(\mathbf{x}\right)={\displaystyle \frac{6\pi \mu R_p^2}{h}}+O\left(1\right)$$

(90)

for small gaps h between sphere and planar wall. At larger distances, Lorentz [60] found

$${\eta }_{3,3}\left(\mathbf{x}\right)=6\pi \mu R_p\left(1+{\displaystyle \frac{9}{8}}{\displaystyle \frac{R_p}{h}}\right)$$

(91)

and, since $9/8\approx 1$, a reasonable approximation over the whole range of h values is given by

$${\eta }_{3,3}\left(\mathbf{x}\right)=6\pi \mu R_p\left(1+{\displaystyle \frac{R_p}{h}}\right)$$

(92)

The singularity of ${\eta }_{3,3}\left(\mathbf{x}\right)$ at the wall implies that a particle, initially at $h=h_0$ and moving towards the solid surface due to the action of a constant force $F_g$, say gravity, reaches the surface in an infinity time $t_{\mathrm{surf}}$, being the integral

$$t_{\mathrm{surf}}={\int }_{h_0}^0{\displaystyle \frac{{\eta }_{3,3}\left(h\right)}{F_g}}dh\to \infty$$

(93)

divergent to infinity. The characteristic time $t_{\mathrm{surf}}$ is referred to as the wall touching time. This result poses severe problems in predicting the kinetics of coalescence and deposition of dispersed particles from hydrodynamic theories based on the Stokes model Equation (5). To complete the picture, let us consider the scaling behavior of the remaining resistance coefficients. As regards the other coefficients, they are all logarithmically singular as $O(log(h\left)\right)$ with the exception of the rotational resistance around the normal axis to the surfaces attaining a finite value close to the wall, i.e., ${\eta }_{6,6}\left(h\right)=O\left(1\right)$ as $h\to 0$. Singularities at a contact point between surfaces are typical of Stoke’s flows due to the no-slip conditions at solid boundaries. For example, in contact line motions [61] and for flows near a corner [62] it has be found that unphysical singularities can be eliminated or mollified by introducing slippage at solid boundaries [63,64]. Before, addressing the effect of slip boundary conditions, it is useful to investigate further the pathologies arising from the classical no-slip hydrodynamic description in dealing with surface phenomena.

Surface Phenomena and Hydrodynamic Singularities

All the surface phenomena (particle coalescence, surface aggregation, adsorption, surface chemical reactions) depending on a transfer mechanism of molecules and particles from the fluid phase onto the surface are deeply influenced by the spatial dependence of the entries of the hydrodynamic resistance matrix. Paradoxes arise due to the singularities of these entries at the solid walls (for incompressible flows, assuming no-slip boundary conditions at the solid boundaries).

In order to highlight these phenomena, it is sufficient to consider a simple problem, namely the pure diffusional motion of solute particles (nanoparticles) in the neighborhood of a solid wall (located at $x_3=x=0$), undergoing at the solid wall a surface chemical reaction characterized by a linear, first-order kinetics. To simplify the analysis, let us assume that far away from the wall, say at $x=L$, the particle concentration is kept fixed, and equal to $c_0$. The problem is thus specified by the parabolic diffusion equation for the particle concentration $c(x,t)$,

$$\frac{\partial c(x,t)}{\partial t}=\frac{\partial }{\partial x}\left(D\left(x\right)\frac{\partial c(x,t)}{\partial x}\right)$$

(94)

equipped with the boundary conditions

$${\left.D\left(x\right)\frac{\partial c(x,t)}{\partial x}=kc(x,t)\right|}_{x=0},c(L,t)=c_0$$

(95)

Consider the steady-state solution $c^{*}\left(x\right)$ of Equation (94). From Equation (94) it follows that,

$$\frac{d{c}^{*}\left(x\right)}{dx}=\frac{B}{D\left(x\right)}$$

(96)

where B is a constant, and thus

$$c^{*}\left(x\right)=A+B{\int }_0^x\frac{d\xi }{D\left(\xi \right)}$$

(97)

where A is a second constant to be determined from the boundary conditions. From Equation (97) it follows that a solution exists provided that $1/D\left(x\right)$ is locally integrable near $x=0$. But this is not the case when the hydrodynamic modeling deriving from incompressible Stokes equations imposing no-slip boundary conditions at the solid walls, as $1/D\left(x\right)\sim 1/x$.

On the other hand, in the case $1/D\left(x\right)$ would be locally integrable near $x=0$, the steady-state solution $c^{*}\left(x\right)$ would attain the expression

$$c^{*}\left(x\right)=c_0\frac{1+k\Phi \left(x\right)}{1+k\Phi \left(L\right)},\Phi \left(x\right)={\int }_0^x\frac{d\xi }{D\left(\xi \right)}$$

(98)

It is clear from the above problem that Stokes’ hydrodynamics applied to transport phenomena coupled to any form of superficial chemical physical processes determines unphysical paradoxes. In the present case, the paradox can be resolved within a classical hydrodynamic formulation by deriving, from more general hydrodynamic conditions, a diffusion coefficient $D\left(x\right)$ that is integrable near $x=0$, i.e., such that $D\left(x\right)\le C/x^{\alpha }$, with $C>0$ and $\alpha <1$. Indeed, this property can be recovered by relaxing the no-slip boundary conditions as addressed in the next Section.

5. Effect of Slip Boundary Conditions

Although the no-slip assumption, $\mathbf{u}=0$, at the interface between a Newtonian fluid and a solid boundary is largely accepted due to its capability of predicting mechanical and hydrodynamical properties involving macroscopic bodies and large-scale systems, the nature of the proper boundary conditions for the tangential velocity at a solid interface has been a long-debated issue over more than two centuries of hydrodynamic research [65]. As an alternative to the no-slip hypothesis, Navier [66] proposed that the tangential stresses at any point on the solid surface should be the same as the stresses at a neighboring internal point of the fluid, providing the boundary condition

$$\beta (\mathbf{u}-\mathbf{v})·(\mathbf{I}-\mathbf{n}\otimes \mathbf{n})=-\mathbf{n}·\mathsf{\tau }·(\mathbf{I}-\mathbf{n}\otimes \mathbf{n})$$

(99)

where $\mathbf{I}$ is the identity matrix, $\mathbf{n}$, the unit normal vector to the surface of the solid, $\beta$ a friction constant and “⊗” indicates dyadic composition, ${(\mathbf{n}\otimes \mathbf{n})}_{i,j}=n_i{n}_j$. Since the isotropic pressure contribution entering $\mathsf{\tau }$ vanishes at the r.h.s of Equation (99), the Navier’s boundary condition (99) can be written as

$$(\mathbf{u}-\mathbf{v})·(\mathbf{I}-\mathbf{n}\otimes \mathbf{n})=\lambda \mathbf{n}·(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T)·(\mathbf{I}-\mathbf{n}\otimes \mathbf{n})$$

(100)

where $\lambda =\mu /\beta$, having the dimension of a length, is the so called slip length. It is evident that the slip length represents a new parameter in modeling Stokes flow. In fact, while the no-slip model describes an idealized fluid/solid interface, the slip model introduces an additional parameter related to the chemical physical interactions at the solid-liquid interface.

Slip phenomena can be distinguished in two main classes: intrinsic slip and apparent slip [67]. The first one is due exclusively to the molecular dynamics at the solid-liquid interface. The second one is due to artifacts at the surfaces that can increase or decrease the slippage, such as the presence of gas bubbles or the surface roughness. However, these artifacts have necessarily a characteristic length below which no-slip conditions applies. In addition, as shown by Cox [58], the nature of the singularity does not depend on the curvature of the approaching surfaces. In considering very small gaps between the surfaces ($h\to 0$), the apparent slip vanishes. Therefore, in the remainder, we will refer solely to the intrinsic (molecular) slip. The occurrence of slippage depends in principle, either on hydrodynamic conditions or on the gap length. In fact, as shown by molecular dynamic simulations [68], the slip length for the system (PA-6,6 oligomer)-graphene increases by increasing the shear rate and by reducing the gap between the solid walls. Being such behavior, due to the molecular rearrangement of the fluid, it is appreciable solely when the gap is small enough to be comparable with the characteristic size of the molecular structure of the fluid. However, at this lengthscale scale, it is not possible to mark a clear distinction between the solid surface and the domain of the fluid used in continuous hydrodynamic models, due to solvation and diffusive phenomena. Therefore, a constant slip length, used in hydrodynamic theories (and in the remainder), should be considered as a parameter emerging from the complexity of the molecular system. A review of typical slip lengths for several fluid-surface systems obtained experimentally is given in [67]. Non vanishing values of the slip length are in the range $\lambda =1÷100$ nm, whereas typical characteristic lengthscales of colloids are $R_p=10÷{10}^4$ nm. Correspondingly, the dimensionless slip length for a colloid attains values in the range $\widehat{\lambda }={10}^{-4}÷10$.

Owing to the linearity of the slip boundary conditions Equation (100), a linear relation between resistance and particle velocity Equation (7) still holds. In [33], a purely mechanical proof of the symmetry of the resistance matrix has been provided, alternative to the thermodynamic approach followed in [69]. However, this proof, as opposed to the thermodynamic one, considers solely no-slip conditions at the solid surfaces. To fill this gap, Appendix A provides a more general proof of the symmetry of the resistance matrix also in the presence of slippage enforcing mechanical arguments.

Either the problems of translating or rotating sphere in an unbounded fluid have been solved by Basset [70] in the presence of the Navier’s slip obtaining the following expressions for the force and the torque acting on the sphere

$${\mathbf{F}}_{\mathrm{hydro}}=-6\pi \mu R_p\left({\displaystyle \frac{1+2\widehat{\lambda }}{1+3\widehat{\lambda }}}\right)\mathbf{v}$$

(101)

$${\mathbf{T}}_{\mathrm{hydro}}=-8\pi \mu R_p^3\left({\displaystyle \frac{1}{1+3\widehat{\lambda }}}\right)\mathsf{\omega }$$

(102)

Being $\widehat{\lambda }={10}^{-4}÷10$, larger spherical particles experience the force and the torque corresponding to no-slip boundary conditions ${\mathbf{F}}_{\mathrm{hydro}}\approx -6\pi \mu R_p\mathbf{v}$ and ${\mathbf{T}}_{\mathrm{hydro}}\approx -8\pi \mu R_p^3\mathsf{\omega }$, whereas for smaller colloids the Basset’s correction becomes necessary.

For a particle in a confined geometry, the slip boundary conditions can be considered either at the surface of the particle ($S_p$), with a slip length ${\lambda }_p$, or at the surface of the walls ($S_w$) with a slip length ${\lambda }_w$. Therefore, the boundary conditions for the flow become

$$\left\{\begin{array}{c}(\mathbf{u}-\mathbf{v})·(\mathbf{I}-\mathbf{n}\otimes \mathbf{n})={\lambda }_p\mathbf{n}·(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T)·(\mathbf{I}-\mathbf{n}\otimes \mathbf{n}),\mathbf{x}\in S_p\hfill \\ \\ (\mathbf{u}-\mathbf{v})·(\mathbf{I}-\mathbf{n}\otimes \mathbf{n})={\lambda }_w\mathbf{n}·(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T)·(\mathbf{I}-\mathbf{n}\otimes \mathbf{n}),\mathbf{x}\in S_w\hfill \end{array}\right.$$

(103)

Spherical Particle Moving Perpendicular to a Planar Wall

Let us to consider the problem of a spherical particle near a planar wall as depicted in Figure 1. Hocking has shown [71] that, in the presence of the same slip on both the surfaces (${\lambda }_p={\lambda }_w=\lambda$), the touching time in Equation (93) becomes finite although the singularity of the transversal resistance at the wall still remains [71]. The functional dependence of the resistance to perpendicular translations on h in the limit $h\to 0$, derived by Hocking using lubrication methods, attains the form

$${\displaystyle \frac{{\eta }_{3,3}(\widehat{h},\widehat{\lambda })}{{\eta }_{\infty }}}={\displaystyle \frac{1}{3\widehat{\lambda }}}\left[\left(1+{\displaystyle \frac{\widehat{h}}{6\widehat{\lambda }}}\right)log\left(1+{\displaystyle \frac{6\widehat{\lambda }}{\widehat{h}}}\right)-1\right]+O\left(1\right)$$

(104)

where $\widehat{h}=h/R_p$. Consequently, the drag force is logarithmically singular at the wall, Equation (91) becomes integrable and the touching time attains a finite value.

A seminalytical expression for the drag force over a sphere translating perpendicularly to a plane has been obtained by Goren [72] over the entire range of positions and for any values of ${\lambda }_w$ and ${\lambda }_p$ using a bispherical coordinate system. However, the author [72] has provided only few numerical values for ${\widehat{\eta }}_{3,3}(\widehat{h},\widehat{\lambda })$ corresponding to relatively large gaps and slip lengths, and considering solely the case ${\lambda }_w={\lambda }_p$. In point of fact, in order to obtain the numeric values for the transversal resistance from the Goren’s solution is in principle necessary to solve an infinite dimensional linear system of equations for any position of the sphere. On the other hand, the infinite linear system, truncated to a finite number N of equations, provides approximate but accurate values for the force experienced by the particle with an error that decreases as N increases, while the number of the equations from a fixed accuracy increases as $\widehat{h}\to 0$ and the slip lengths increase.

In order to obtain the values of ${\eta }_{3,3}(\widehat{h},{\widehat{\lambda }}_p,{\widehat{\lambda }}_w)$ for different slip lengths and at any distance from the wall, we have either performed FEM simulations, the detail of which are reported in the Appendix B, or solved the Goren’s equations up to $N=500$. The data depicted in Figure 2 show that the Hocking asymptotic Equation (104) matches accurately the FEM simulations and the Goren solution in the range of $\widehat{\lambda }=0÷{10}^{-1}$ and that the predicted logarithmic scaling starts to appear for gaps smaller then the slip length ($\widehat{h}\lesssim \widehat{\lambda }$), after a transition zone, where ${\eta }_{3,3}\sim 1/\widehat{h}$.

On the other hand, as depicted in Figure 3, if the no-slip boundary condition is imposed on solely one of the surfaces, the singular scaling ${\eta }_{3,3}\sim 1/\widehat{h}$ remains no matter the value of the slip length imposed on the other surface. In the latter conditions, we can distinguish three different regimes: (i) the scaling ${\eta }_{3,3}\sim 1/\widehat{h}$ for $\widehat{\lambda }\lesssim \widehat{h}\lesssim {10}^{-1}$, (ii) an apparent logarithmic behavior for ${10}^{-2}\widehat{\lambda }\lesssim \widehat{h}\lesssim \widehat{\lambda }$ and (iii) the asymptotic regime where ${\eta }_{3,3}\sim 1/\widehat{h}$ for $\widehat{h}\lesssim {10}^{-2}\widehat{\lambda }$.

To complete the analysis the data in Figure 4 show that, keeping fixed the slip length at one of the surfaces (${\widehat{\lambda }}_i={10}^{-3}$, $i=w,p$) and increasing the slip length on the other ${\widehat{\lambda }}_j$, $j=p,w$, the logarithmic scaling occurs for $\widehat{h}\lesssim {\widehat{\lambda }}_j$. This means that an arbirarily small slip on both the surfaces is sufficient to determine an asymptotic logarithmic scaling of the transversal resistance and thus the occurrence of a finite value of the touching time.

6. Fluid Inertial Effects

So far we have considered fluid-particle interaction within the instantaneous Stokes regime, neglecting fluid inertia. In point of fact, while the Reynolds number is smaller that 1 in most of the microfluidic applications, this is not the case of the product of the Reynolds number times the Strouhal number, which is order of 1 or higher due to the high frequency of the thermal fluctuations. This means that a more accurate description of particle motion at short time and length scales would involve a hydrodynamic description of fluid inertia, which in the present case corresponds to the time-dependent Stokes regime [69],

$$\rho \frac{\partial \mathbf{u}}{\partial t}=\mu {\nabla }^2\mathbf{u}-\nabla p$$

(105)

with $\nabla ·\mathbf{u}=0$.

It is well known that the force ${\mathbf{F}}_{\mathrm{hydro}}\left(t\right)$ exerted by a fluid with density $\rho$ and viscosity $\mu$ on a spherical particle of radius $R_p$ moving with velocity $\mathbf{v}$ in a still fluid, can be expressed in the Laplace domain ${\widehat{\mathbf{F}}}_{\mathrm{hydro}}\left(s\right)=L\left[{\mathbf{F}}_{\mathrm{hydro}}\left(t\right)\right]={\int }_0^{\infty }{e}^{-st}{\mathbf{F}}_{\mathrm{hydro}}\left(t\right)dt$ as [34]

$${\widehat{\mathbf{F}}}_{\mathrm{hydro}}\left(s\right)=-6\pi \mu R_p\widehat{\mathbf{v}}\left(s\right)-6\pi \sqrt{\mu \rho }{R}_p^2\frac{1}{\sqrt{s}}\left(s\widehat{\mathbf{v}}\left(s\right)\right)-\frac{2}{3}\pi R_p^3\rho \left(s\widehat{\mathbf{v}}\left(s\right)\right)$$

(106)

The first term at the r.h.s of Equation (106) is the Stokesian friction factor, the second one corresponds in the time domain to the convolutional Basset force

$${\mathbf{F}}_{\mathrm{Basset}}\left(t\right)=-6\sqrt{\pi \mu \rho }{R}_p^2{\int }_0^t\frac{1}{\sqrt{t-\tau }}\left(\frac{d\mathbf{v}\left(\tau \right)}{d\tau }+\mathbf{v}\left(0\right)\delta \left(\tau \right)\right)d\tau$$

(107)

and the third term is the added-mass contribution

$${\mathbf{F}}_{\mathrm{am}}\left(t\right)=-m_a\left(\frac{d\mathbf{v}\left(t\right)}{dt}+\mathbf{v}\left(0\right)\delta \left(t\right)\right),m_a=\frac{V_p\rho }{2}$$

(108)

where $V_p$ is the particle volume, equal to half the mass of the fluid occupying the particle volume. Physically, the added mass contribution corresponds to the back action on the particle of the correlated motion of nearby fluid elements originating from the particle movement within the fluid [73]. In confined geometries, the added mass ${\mathbf{m}}_a\left(\mathbf{x}\right)$ becomes a tensorial quantity dependent on particle position. For instance, for a spherical particle near a solid planar wall,

$${\mathbf{m}}_a\left(\mathbf{x}\right)=\left(\begin{array}{ccc}{m}_1\left(h\right)& 0& 0\\ 0& m_1\left(h\right)& 0\\ 0& 0& m_2\left(h\right)\end{array}\right)$$

(109)

where $m_1\left(h\right)$, and $m_2\left(h\right)$, corresponding to the parallel and transversal added masses, are smooth functions of the particle distance from the wall, as depicted in Figure 5 attaining a finite value for $h=0$.

Similarly, the Basset contribution attains in confined geometries a tensorial, position-dependent character,

$${\mathbf{F}}_{\mathrm{Basset}}={\int }_0^t\mathbf{B}(t-\tau ,\mathbf{x})\frac{d\mathbf{v}\left(\tau \right)}{d\tau }d\tau$$

(110)

where $\mathbf{B}(t,\mathbf{x})={\left(B_{i,j}(t,\mathbf{x})\right)}_{i,j=1}^3$ and $\mathbf{x}=\mathbf{x}\left(t\right)$. There are very few works on the characterization of the Basset force in confined geometries [74], and the detailed study of fluid inertial effects in microchannels and in the presence of solid walls represents an almost virgin field of theoretical and experimental investigation.

Gathering these contributions, and considering also the action of a potential $U\left(\mathbf{x}\right)$ and of a flow force ${\mathbf{F}}_{\mathrm{flow}}=\eta \left(\mathbf{x}\right){\mathbf{u}}^{\left(p\right)}\left(\mathbf{x}\right)$, deriving e.g., by a (pressure-drive) flow in the fluid, the equations of motion for a microparticle become,

$$\begin{array}{cc}\hfill \frac{d\mathbf{x}}{dt}& =\mathbf{v}\\ \hfill (m+{\mathbf{m}}_a\left(\mathbf{x}\right))\frac{d\mathbf{v}}{dt}& =& -\mathsf{\eta }\left(\mathbf{x}\right)(\mathbf{v}-{\mathbf{u}}^{\left(p\right)}\left(\mathbf{x}\right))-\nabla U\left(\mathbf{x}\right)\hfill \\ & -& {\int }_0^t\mathbf{B}(t-\tau ,\mathbf{x})\left(\frac{d\mathbf{v}\left(\tau \right)}{d\tau }-\frac{d{\mathbf{u}}^{\left(p\right)}\left(\mathbf{x}\left(\tau \right)\right)}{d\tau }\right)+{\mathbf{F}}_{\mathrm{stocha}}(\mathbf{x},t)\hfill \end{array}$$

(111)

where ${\mathbf{F}}_{\mathrm{stocha}}(\mathbf{x},t)$ is the stochastic fluctuational contributions, to be determined by enforcing fluctuation-dissipation relations. But also this aspect, owing to the explicit dependence of ${\mathbf{F}}_{\mathrm{stocha}}(\mathbf{x},t)$ on the position $\mathbf{x}$, represents a challenging problem in statistical physics, especially if one is interested in determining the velocity autocorrelation function. Fortunately, for $t\gg t_{\mathrm{diss}}$, corresponding to typical conditions in microfluidic applications, the fluid inertia can be neglected, and the stochastic description of particle motion can be based on the theory addressed in Section 2 and Section 3.

7. Concluding Remarks

Microfluidics is not only an emerging field as regards practical applications but it provides a fertile playground for analyzing fundamental physical questions. This is because particle motion at microscale is controlled by thermal stochastic fluctuations, the understanding of which would not only improve the performances of micofluidic devices (e.g. for separation purposes), but also could test the validity of classical physical theories (hydrodynamics) down to the scale where quantum effects should start to appear (going down to the nanoscale).

The fundamental peculiarity of microhydrodynamics in confined geometries is that fluid-particle interactions depend on the position and, owing to fluctuation-dissipation relations, these spatial nonuniformities are transferred to mass-transport properties (diffusion) giving rise to delicate theoretical issues and paradoxes.

As long as fluid inertia is negligible, it is still possible to derive a rather complete characterization of thermal fluctuations, and to perform a simplification of the equations of motion, eliminating the fast velocity variables, thus reducing the statistical description to the Fokker-Planck equation for the particle spatial density corresponding to a classical advection-diffusion equation in physical space and time. However, even in this case, transport paradoxes arise due to the singularity of the hydrodynamic resistances in the neighbourhood of a solid wall.

We have thoroughly analyzed how the lack of integrability of mass-transport models (diffusion equation) in the presence of a surface chemical reaction depend on the simplifying assumption of no-slip boundary conditions. The inclusion of slippage effects at all the solid surfaces (walls and particle) transforms the non-integrable $1/h$-singularity in the transversal resistance coefficient ${\eta }_{3,3}\left(h\right)$ near a planar wall into a logarithmic singularity (∼$-logh$) resolving the above mentoned integrability problem. Hovewer, there are physical reasons to conjecture that a more refined modeling of fluid-particle interactions (including compressibility and acoustic effects) could completely eliminate the singularity in ${\eta }_{3,3}\left(h\right)$ at the wall. This issue will be addressed in forthcoming works.

The inclusion of fluid inertia in confined microfluidics opens up other interesting and fundamental problems. The theoretical and numerical estimate of the Basset force even in simple confined systems is a basic ingredient in order to tackle the formulation of fluctuation-dissipation relations, and the efficient determination of the fluctuational force ${\mathbf{F}}_{\mathrm{stocha}}(\mathbf{x},t)$, similarly to what developed in the free space, in order to determine the short-time behavior of the velocity autocorrelation functions. In this framework, the development of more elaborate statistical physical methods should proceed hand-in-hand with the accurate experimental analysis of particle motion at short time and length scales near a solid surface, following and extending the experimental investigations recently performed for Brownian particles in the free space [75,76].