1 Introduction

Hydrodynamics involves the solution of a vectorial evolution equation, such as the parabolic Navier–Stokes equations for Newtonian fluids. Two main sources of difficulties are involved with the mathematical solution of hydrodynamic models: its nonlinear nature, owing to fluid inertia, and its vectorial character, as the velocity is a vector field [1]. Nonlinearities arise in many scalar transport models, such as reaction-diffusion equation in a chemical mixture involving N different chemical species: if \(\mathbf{c}(\mathbf{x},t)=(c_1(\mathbf{x},t),\dots ,c_N(\mathbf{x},t))\) is the vector-valued function of their (molar) concentrations, the evolution of this chemical system in a quiescient liquid is described by the parabolic model

$$\begin{aligned} \frac{\partial \mathbf{c}(\mathbf{x},t)}{\partial t}= {{\mathcal {D}}} \, \nabla ^2 \mathbf{c}(\mathbf{x},t)+ \mathbf{r}(\mathbf{c}(\mathbf{x},t)) \end{aligned}$$
(1)

where \({{\mathcal {D}}}=\text{ diag }(D_{1},\dots ,D_N)\) is the diffusivity tensor admitting diagonal form and constant entries (for simplicity of discussion), and \(\mathbf{r}(\mathbf{c})\) is the contribution in the molar balance of the chemical reactions occurring in the mixture. A mathematical analysis of these models has been extensively developed, as concerns their bifurcational properties, at least in some simple but paradigmatic cases [2].

The geometric vectorial nature of the Navier–Stokes equations copes with nonlinearities adding a further level of complexity to hydrodynamic problems. In point of fact, still limiting the analysis to flows characterized by small, and even vanishingly small Reynolds numbers,—occurring in the majority of microfluidic applications or in the analysis of Brownian motion of colloidal particles—so that the Navier–Stokes equations can be conveniently approximated with reliable accuracy by means of the unsteady Stokes model

$$\begin{aligned} \rho \, \frac{\partial \mathbf{v}(\mathbf{x},t)}{\partial t}= \mu \,\nabla ^2 \mathbf{v}(\mathbf{x},t) - \nabla p(\mathbf{x},t) \end{aligned}$$
(2)

where \(\rho\) and \(\mu\) are respectively the fluid density and viscosity, and \(\mathbf{v}(\mathbf{x},t)\), \(p(\mathbf{x},t)\) the velocity and pressure field, the mere assumption of incompressibility,

$$\begin{aligned} \nabla \cdot \mathbf{v}(\mathbf{x},t) =0 \end{aligned}$$
(3)

changes completely the qualitative properties of the parabolic hydrodynamic model (2), determining the occurrence of novel (and sometimes unexpected) features that find no counterpart in the scalar, albeit vector-valued, model Eq. (1).

This is because, the vector Laplacian operator entering Eq. (2) does possess geometric properties that are absent for the Laplacian operator acting on scalar functions as in Eq. (1), which are related to the decomposition

$$\begin{aligned} \nabla ^2 \mathbf{v}(\mathbf{x},t) = - \nabla \times \left[ \nabla \times \mathbf{v}(\mathbf{x},t) \right] + \nabla \left[ \nabla \cdot \mathbf{v}(\mathbf{x},t) \right] \end{aligned}$$
(4)

thus resulting in the double and repeated action of the curl operator with reversed sign when the incompressibility condition (3) is enforced. The splitting of the Laplacian in two operators Eq. (4) is in one-to-one correspondence with the Helmholtz decomposition of the space of square summable vector fields in \({{\mathbb {R}}}^3\) into two orthogonal subspaces possessing respectively vanishing divergence and vanishing curl. This issue and its physical implications are addressed in paragraph 4.1 in connection with the Helmholtz decomposition of vectorial impulsive functions.

This observation can be reformulated as follows: the vectorial nature of the Navier–Stokes equations, once coupled with the constraint associated with the solenoi-dal nature of the velocity field, determines long-distance spatial correlations that have no counterpart in any parabolic scalar model such as Eq. (1). We refer to Sect. 4 for a thorough analysis of this issue.

Two milestones have marked the birth of the hydrodynamic theory of low Reynolds number flows: the introduction and explicit derivation of the hydrodynamic Green function for the Stokes flow due to Oseen in the three-dimensional free propagation (i.e., the derivation of the Oseen tensor) [3], and the mathematical formulation of the hydrodynamic theory due to Olga Ladyzhenskaya, which highlighted in its seminal contribution to the field [4] the mathematical consequences of the vectorial nature of the Navier–Stokes equations, the relationship with the Helmholtz decomposition of a vector field, and the necessity of a pressure gauge in order to properly define the spectral properties of the vector Laplacian in incompressible flows (i.e., the so called Ladyzhenskaya theorems).

The derivation of the Oseen tensor and the subsequent analysis on the singularity representation of Stokes flows paved the way towards a robust theory of low-Reynolds number hydrodynamics, leading naturally to the formulation of efficient boundary element methods for the numerical solution of low-Reynolds number problems [5,6,7].

The aim of the present article is to provide a simple method for deriving the hydrodynamic Green function for unsteady incompressible Stokes flows—which of course can be applied under steady Stokes conditions—and derive its physical consequences. The method is essentially based on the resolution of two scalar problems, although one of which involves vector-valued functions. The leading idea is very simple as it involves the free propagation of the velocity field without any vectorial condition, as Eq. (2), on its solenoidal nature, restoring this property upon a gradient gauge. This way of approaching the hydrodynamic problem permits to highlight the nonlocality of the incompressibility constraint (3) in unsteady Stokes regime, deriving from the solution of a Poisson equation in which time enters solely as a parameter. The application of this method to the unsteady Stokes equation is particularly interesting, not only because it provides a simple and closed-form expression for the Green function, which has been already studied in several articles [6, 8,9,10], but mainly because the application of the gradient gauge, characterizing this formulation, permits to understand the physical paradox associated with the large-scale power-law scaling of the Green function with the distance from the point source at any time \(t>0\), and to relate it to the infinite velocity of propagation characterizing any incompressible hydrodynamic models.

The article is organized as follows. Section 2 introduces the method using the derivation of the Oseen tensor as a test case. The same approach is applied to the unsteady Stokes equation in Sect. 3, and a simple closed-form expression for the Green function obtained. Section 4 develops the qualitative analysis of the functional form of the hydrodynamic Green function in unsteady Stokes conditions, its regularity structure at infinity and the role of incompressibility on the unbounded propagation velocity of stresses (Stokesian propagation paradox). The propagation paradox associated with incompressibility is also explained directly from the nature of the forcing term defining the hydrodynamic Green functions, by considering the Helmholtz decomposition of an impulsive vector-valued Dirac’s delta function. Section 5 develops the boundary integral formulation deriving from the gradient gauge approach for hydrodynamic problems in the presence of solid boundaries.

2 Explanation of the method: calculation of the Oseen tensor

Although it is possible to derive the Oseen tensor for Stokes flows in several different ways [6, 7, 12, 13], at least four as observed by Lisicki [14], we use this prototypical problem to set up a simple and general approach for deriving hydrodynamic Green functions which can be easily extended to more complex and interesting problems.

Consider the Green function for the Stokes problem in an incompressible flow

$$\begin{aligned} \mu \, \nabla ^2 \mathbf{v}(\mathbf{x}) - \nabla p(\mathbf{x}) + \mathbf{f}_0 \, \delta (\mathbf{x}-\mathbf{x}_0) =0 \end{aligned}$$
(5)
$$\begin{aligned} \nabla \cdot \mathbf{v}(\mathbf{x}) =0 \end{aligned}$$
(6)

\(\mathbf{x} \in {{\mathbb {R}}}^3\), where \(\mathbf{f}_0\) is a vector-valued constant. This set of equations is equipped with the regularity condition at infinity, \(\mathbf{v} \rightarrow 0\) for \(|\mathbf{x}|\rightarrow \infty\), and with some regularity constraints in \({{\mathbb {R}}}^3\). From summability conditions, we have to enforce that \(\mathbf{v}(\mathbf{x})\) should diverge in the neighbourhood of any point \(\mathbf{x}^\prime \in {{\mathbb {R}}}^3\) slower than \(1/|\mathbf{x}-\mathbf{x}^\prime |^{\alpha }\) with \(\alpha <2\) [15].

The main idea is to decompose the velocity \(\mathbf{v}(\mathbf{x})\) into two contributions

$$\begin{aligned} \mathbf{v}(\mathbf{x}) = \mathbf{v}^\prime (\mathbf{x}) + \nabla \phi (\mathbf{x}) \end{aligned}$$
(7)

The vector \(\mathbf{v}^\prime (\mathbf{x})\) accounts for the viscous dynamics, i.e., it is the solution of the equation

$$\begin{aligned} \mu \, \nabla ^2 \mathbf{v}^\prime (\mathbf{x}) + \mathbf{f}_0 \, \delta (\mathbf{x}-\mathbf{x}_0) =0 \end{aligned}$$
(8)

without imposing any further condition on its divergence. From what discussed in the Introduction, Eq. (8) corresponds to a scalar Poisson equation applied to the vector-valued quantity \(\mathbf{v}^\prime (\mathbf{x})\). The velocity field \(\mathbf{v}^\prime (\mathbf{x},t)\) in itself does not have a direct hydrodynamic meaning.

Conversely, the gradient gauge \(\nabla \phi (\mathbf{x})\) takes care of Eq. (6), implying, from Eq. (7), that

$$\begin{aligned} \nabla ^2 \phi (\mathbf{x}) = - \nabla \cdot \mathbf{v}^\prime (\mathbf{x}) \end{aligned}$$
(9)

The hydrodynamic problem is completely solved by determining \(\mathbf{v}^\prime (\mathbf{x})\) and \(\nabla \phi (\mathbf{x})\) from Eqs. (8–9), and the value of the pressure \(p(\mathbf{x})\) follows from Eqs. (5) and (7) as

$$\begin{aligned} p(\mathbf{x}) = p_0 + \mu \, \nabla ^2 \phi (\mathbf{x}) \end{aligned}$$
(10)

where \(p_0\) is any constant value.

Consider Eq. (8). This is a Poisson equation in the free space \({{\mathbb {R}}}^3\), thus admitting \(E(\mathbf{x})= - 1/4 \, \pi |\mathbf{x}|\) as its fundamental solution, \(\nabla ^2E(\mathbf{x})=\delta (\mathbf{x})\), (Green function). Consequently, the solution of Eq. (8) takes the form

$$\begin{aligned} \mathbf{v}^\prime (\mathbf{x})= & {} \frac{\mathbf{f}_0}{4 \, \pi \, \mu } \int _{{\mathbb R^3}} \frac{\delta (\varvec{\xi }-\mathbf{x}_0)}{|\mathbf{x}-\varvec{\xi }|} \, d\varvec{\xi }+ \mathbf{v}_\mathrm{hom}^\prime (\mathbf{x}) \nonumber \\= & {} \frac{1}{4 \, \pi \, \mu } \frac{\mathbf{f}_0}{|\mathbf{x}-\mathbf{x}_0|} + \mathbf{v}_\mathrm{hom}^\prime (\mathbf{x}) \end{aligned}$$
(11)

where \(\mathbf{v}_\mathrm{hom}^\prime (\mathbf{x})\) is any solution of the Laplace equation in the free space. These homogeneous solutions of the Laplace equation can be written as the combination of the elementary solutions \(\mathbf{v}_\mathrm{hom}^\prime (\mathbf{x}) = A \, \nabla (|\mathbf{x} -\mathbf{x}^\prime |^{-1})\) for generic \(\mathbf{x}^\prime\), where A is an arbitrary constant. Specifically, for \(\mathbf{x}^\prime =\mathbf{x}_0\), this homogeneous solution provides the velocity source singularity \(\nabla (|\mathbf{x} -\mathbf{x}_0|^{-1})\). However, these contributions do not possess the required regularity properties near \(\mathbf{x}^\prime\) (as the velocity diverges as \(|\mathbf{x}-\mathbf{x}^\prime |^{-2})\), and consequently they should be discarded, i.e., the multiplicative constant A should be vanishing.

The calculation of the divergence of \(\mathbf{v}^\prime (\mathbf{x})\) provides

$$\begin{aligned} \nabla \cdot \mathbf{v}^\prime (\mathbf{x})= \frac{1}{4 \, \pi \, \mu } \nabla _x \left( \frac{1}{|\mathbf{x}-\mathbf{x}_0|} \right) \cdot \mathbf{f}_0 \end{aligned}$$
(12)

Equation (9) is again a free-space Poisson equation, the solution of which, making use of Eq. (12) reads

$$\begin{aligned} \phi (\mathbf{x})= - \frac{1}{16 \, \pi ^2 \, \mu } \int _{{{\mathbb {R}}}^3} \frac{1}{|\mathbf{x}-\varvec{\xi }|} \, \nabla _\xi \left( \frac{1}{|\varvec{\xi } - \mathbf{x}_0|} \right) \cdot \mathbf{f}_0 \, d \varvec{\xi } \end{aligned}$$
(13)

Equations (11) and (13) formally solve the mathematical problem of estimating the Stokes Green function. Nevertheless, it is convenient to explicit further the formal solution for \(\phi (\mathbf{x})\). Set \(\mathbf{r}=\mathbf{x}-\mathbf{x}_0\). The scalar gauge \(\phi (\mathbf{r})\) can be written as

$$\begin{aligned} \phi (\mathbf{r}) = \frac{1}{4 \, \pi \, \mu } \, u(\mathbf{r}) \end{aligned}$$
(14)

where \(u(\mathbf{r})\) is the solution of the Poisson equation

$$\begin{aligned} \nabla ^2 u(\mathbf{r}) = \frac{\mathbf{r} \cdot \mathbf{f}_0}{r^3} \end{aligned}$$
(15)

where \(r= |\mathbf{r}|\), that can be sought in the form

$$\begin{aligned} u(\mathbf{r})= b \, \frac{\mathbf{r} \cdot \mathbf{f}_0}{r^\alpha } \end{aligned}$$
(16)

with constants b and \(\alpha\) to be determined. Substituting the functional form (16) into Eq. (15) and developing elementary operations, one gets \(\alpha =1\), \(b=-1/2\). These calculations involve essentially the properties: \(\nabla r =\mathbf{r}/r\), \(\nabla \cdot \mathbf{r}=3\), \(\nabla \left( \mathbf{r} \cdot \mathbf{f}_0 \right) =\mathbf{f}_0\). Therefore, the scalar gauge \(\phi (\mathbf{r})\) takes the simple form

$$\begin{aligned} \phi (\mathbf{r})=- \frac{1}{8 \, \pi \, \mu } \frac{\left( \mathbf{r} \cdot \mathbf{f}_0 \right) }{r} \end{aligned}$$
(17)

Considering that

$$\begin{aligned} \nabla \left( \frac{(\mathbf{r} \cdot \mathbf{f}_0)}{r} \right) = \frac{\mathbf{f}_0}{r}- \frac{\mathbf{r} \left( \mathbf{r} \cdot \mathbf{f}_0 \right) }{r^3} \end{aligned}$$
(18)

the velocity field \(\mathbf{v}(\mathbf{r})\) is given by

$$\begin{aligned} \mathbf{v}(\mathbf{r})= & {} \frac{1}{4 \, \pi \, \mu } \frac{\mathbf{f}_0}{r} - \frac{1}{8 \, \pi \, \mu } \nabla \left( \frac{(\mathbf{r} \cdot \mathbf{f}_0)}{r} \right) \nonumber \\= & {} \frac{1}{8 \, \pi \, \mu } \left( \frac{\mathbf{f}_0}{r} + \frac{ \mathbf{r} \, ( \mathbf{r} \cdot \mathbf{f}_0)}{r^3} \right) \end{aligned}$$
(19)

Equation (19) permits to defines the entries \(G_{i,j}(\mathbf{r})\) of the tensorial Green function (the Oseen tensor), as

$$\begin{aligned} \mathbf{v}(\mathbf{x}) = {\mathbf {G}}(\mathbf{x}-\mathbf{x}_0) \cdot {\mathbf {f}}_{0}= \sum _{j=1}^3 G_{i,j}(\mathbf{x}-\mathbf{x}_0) \, f_{0,j}, \end{aligned}$$
(20)

with

$$\begin{aligned} G_{i,j}(\mathbf{x}) = \frac{1}{8 \, \pi \, \mu } \left( \frac{\delta _{i,j}}{r}+ \frac{x_i \, x_j}{r^3} \right) \end{aligned}$$
(21)

where \(\mathbf{x}= \sum _i x_i \mathbf{e}_i\), \(\{ \mathbf{e}_i \}_{i=1}^3\) is an orthonormal Cartesian base, \(r=|\mathbf{x}|\), and \(\delta _{i,j}\) are the Kronecker symbols. From Eqs. (10) and (12) the pressure is given by

$$\begin{aligned} p(\mathbf{x}) - p_0 = {\mathbf {g}}(\mathbf{x}-\mathbf{x}_0) \cdot {\mathbf {f}}_0; \mathrm {with} g_i(\mathbf{x}) = \frac{x_i}{4 \, \pi \, r^3} \end{aligned}$$
(22)

where \(p_0\) is any constant.

The above analysis is independent of the coordinate system chosen to represent the spatial variable \(\mathbf{x}\). Nevertheless, Eqs. (5), (8), as well as all the equations of hydrodynamic Green function theory (see e.g. the monographs [5,6,7]) are not written in a fully consistent covariant way, and this may generate confusion whenever a reference system different from a Cartesian one is chosen to represent componentwise the velocity field \(\mathbf{v}\). A fully covariant generalization of hydrodynamic Green function theory can be achieved by applying to it the methods of bitensor analysis originally developed by Synge [16] introducing the concept of world-function in general relativity [17]. This extension will be addressed in a future work.

3 Green function for the unsteady Stokes problem

The approach outlined in the previous section can be applied on equal footing to determine analytically, and in a very simple way, the Green function for the unsteady Stokes problem in the free space. In this case, instead of Eq. (5) one has

$$\begin{aligned} \rho \, \frac{\partial \mathbf{v}(\mathbf{x},t)}{\partial t} = \mu \, \nabla ^2 \mathbf{v}(\mathbf{x},t) - \nabla p(\mathbf{x},t) + \mathbf{f}_0 \, \delta (\mathbf{x}-\mathbf{x}_0) \, \delta (t) \end{aligned}$$
(23)

while incompressibility condition Eq. (6) still holds for the unsteady velocity field \(\mathbf{v}(\mathbf{x},t)\). The initial condition is \(\mathbf{v}(\mathbf{x},0)=0\). Enforcing the decomposition (7), where in the present case all the field \(\mathbf{v}(\mathbf{x},t)\), \(\mathbf{v}^\prime (\mathbf{x},t)\), \(\phi (\mathbf{x},t)\) depends explicitly on time t, the viscous component \(\mathbf{v}^\prime (\mathbf{x},t)\) is the solution of equation

$$\begin{aligned} \rho \, \frac{\partial \mathbf{v}^\prime (\mathbf{x},t)}{\partial t} = \mu \, \nabla ^2 \mathbf{v}^\prime (\mathbf{x},t) + \mathbf{f}_0 \, \delta (\mathbf{x}-\mathbf{x}_0) \, \delta (t) \end{aligned}$$
(24)

with \(\mathbf{v}^\prime (\mathbf{x},0)=0\), while \(\phi (\mathbf{x},t)\) solves the Poisson problem (9) by accounting the explicit dependence on time t of \(\mathbf{v}^\prime (\mathbf{x},t)\). In this elliptic problem time t plays the role of a parameter.

As regards the boundary conditions, regularity at infinity implies that \(\mathbf{v}(\mathbf{x},t)\) should decay to zero for \(|\mathbf{x}| \rightarrow \infty\). A further discussion on the regularity condition at infinity is addressed in the next section.

The solution of Eq. (24) involves the classical scalar Gaussian heat kernel,

$$\begin{aligned} G_D(r,t)= \frac{e^{-r^2/4 \, \nu \, t}}{(4 \, \pi \, \nu \, t)^{3/2}} \end{aligned}$$
(25)

where \(\mathbf{r}=\mathbf{x}-\mathbf{x}_0\), \(r=|\mathbf{r}|\) and \(\nu =\mu /\rho\)

$$\begin{aligned} \mathbf{v}^\prime (\mathbf{x},t)= G_D(r,t) \, \widetilde{\mathbf{f}}_0 = \frac{1}{(4 \, \pi \, \nu \, t)^{3/2}} e^{-r^2/4 \, \nu \, t} \, \widetilde{\mathbf{f}}_0 \end{aligned}$$
(26)

where \(\widetilde{\mathbf{f}}_0 = \mathbf{f}_0/\rho\). Conversely, \(\phi (\mathbf{x},t)\) is the solution of the equation,

$$\begin{aligned} \nabla ^2 \phi (\mathbf{x},t)= & {} - \nabla \cdot \mathbf{v}^\prime (\mathbf{x},t) \nonumber \\= & {} - \nabla \left[ \frac{1}{(4 \, \pi \, \nu \, t)^{3/2}} e^{-r^2/4 \, \nu \, t} \right] \cdot \widetilde{\mathbf{f}}_0 \end{aligned}$$
(27)

Consider a function \(f \left( y \right) = \text{ erf } \left( y \right) /y\), of the radius \(y=|\mathbf{y}|\), \(\text{ erf } \left( y \right) = \frac{2}{\sqrt{\pi }} \int _0^y e^{- \eta ^2} d \eta\) being the error function, where \(\mathbf{y}\) is a normalized position vector to be specified below. It satisfies the property

$$\begin{aligned} \nabla _y^2 \left( \frac{1}{y} \text{ erf } \left( y \right) \right) = - \frac{4}{\sqrt{\pi }} e^{- y^2} \end{aligned}$$
(28)

where the Laplacian \(\nabla _y^2\) refers to the \(\mathbf{y}\)-coordinates. Applying this last equality to Eq. (27), with \({\mathbf {y}} = {\mathbf {r}} / \sqrt{4 \nu t}\), we find:

$$\begin{aligned} \phi (\mathbf{x},t) = \left. \frac{1}{16 \pi \nu t} \nabla _y \left[ \frac{1}{y} \text{ erf } \left( y \right) \right] \cdot \widetilde{\mathbf{f}}_0 \right| _{\mathbf{y}=(\mathbf{x}-\mathbf{x}_0)/\sqrt{4 \nu t}} \end{aligned}$$
(29)

Therefore, we may conclude that the velocity field reads:

$$\begin{aligned} {\mathbf {v}} \left( {\mathbf {x}}, t \right)= & {} {\mathbf {v}}^{\prime } \left( {\mathbf {x}}, t \right) + \nabla \phi ({{\mathbf {x}}},t) = \frac{1}{(4 \pi \nu t)^{3/2}} \left\{ e^{-y^2} {\mathbf {I}} \right. \nonumber \\+ & {} \left. \left. \frac{\sqrt{\pi }}{4} \nabla _y \nabla _y \left[ \frac{1}{y} \text{ erf } \left( y \right) \right] \right\} \cdot \widetilde{\mathbf{f}}_0 \right| _{\mathbf{y}=(\mathbf{x}-\mathbf{x}_0)/\sqrt{4 \nu t}} \end{aligned}$$
(30)

which coincides with Ladyzhenskaya’s solution [4], as it has been expressed in Mauri and Rubinstein [9] (see also [10], and the solutions in [8] and [11]).

This expression for the velocity field can be written in an alternative form applying the identity

$$\begin{aligned} \nabla _y \nabla _y f(y) = \frac{1}{y} \frac{d f(y)}{dy} {\mathbf {I}} + \frac{1}{y} \frac{d}{dy} \left( \frac{1}{y} \frac{d f(y)}{dy} \right) {\mathbf {y}} {\mathbf {y}} \end{aligned}$$
(31)

to \(f(y) = \text{ erf } \left( y \right) /y\) - in Eq. (31), \(\mathbf{I}\) is the identity operator and \(\mathbf{y}{} \mathbf{y}\) the dyadic tensor, the entries of which are \(y_i \, y_j\), often referred as \(\mathbf{y} \otimes \mathbf{y}\) -, obtaining

$$\begin{aligned} \mathbf{v}(\mathbf{x},t)= \frac{1}{(4 \pi \nu t)^{3/2}} \left. \left[ \, a(y) \, {\mathbf {I}} + b(y) \, \frac{\mathbf{y} \mathbf{y}}{y^2} \right] \cdot \widetilde{\mathbf{f}}_0 \right| _{\mathbf{y}=(\mathbf{x}-\mathbf{x}_0)/\sqrt{4 \nu t}} \end{aligned}$$
(32)

where the auxiliary scalar functions a(y) and b(y) are defined as:

$$\begin{aligned} a(y) = e^{-y^2} -h(y) \, , \qquad b(y)=-e^{-y^2} + 3 \, h(y) \end{aligned}$$
(33)

where

$$\begin{aligned} h(y)= \frac{\sqrt{\pi } \, \text{ erf }(y)}{4 y^3}- \frac{1}{2 y^2} e^{-y^2} \end{aligned}$$
(34)

Observe from Eqs. (33–34), the short-time and long-distance scaling \(a(y) \sim b(y) \sim 1/y^3\) corresponding to the terms involving the error function.

As for the pressure field, referred to a constant reference value, as in Eq. (10), from the governing equations we obtain:

$$\begin{aligned} p \left( {\mathbf {x}}, t \right) = \mu \nabla ^2 \phi (\mathbf{x},t) - \rho \frac{\partial \phi (\mathbf{x},t)}{\partial t} \end{aligned}$$
(35)

which generalizes Eq. (10). Considering that from Eq. (27), we have:

$$\begin{aligned} \nabla ^2 \phi (\mathbf{x},t) = \left. - \nabla G_D \left( y, t \right) \cdot \widetilde{\mathbf{f}}_0 \, \right| _{y=r/\sqrt{4 \nu t}} \end{aligned}$$
(36)

from Eq. (35) we obtain:

$$\begin{aligned} \nabla ^2 p = \nabla \left[ \rho \frac{\partial G_D}{\partial t} - \mu \nabla ^2 G_D \right] = \nabla \delta \left( {\mathbf {x}} \right) \delta \left( t \right) \cdot \widetilde{\mathbf{f}}_0. \end{aligned}$$
(37)

This result can also be obtained directly by taking the divergence of Eq. (23). Now, considering that by symmetry the solution must be of the form \(p \left( {\mathbf {x}}, t \right) = {\mathbf {a}} \left( {\mathbf {x}} \right) \cdot \widetilde{\mathbf{f}}_0 \delta \left( t \right)\) with \(\nabla ^2 {\mathbf {a}} \left( \mathbf{x} \right) = \nabla \delta \left( {\mathbf {x}} \right)\), and reminding that the Green function of the Poisson equation in the free space, \(\nabla ^2 E(\mathbf{x}) = \delta \left( \mathbf{x} \right)\), is \(E(\mathbf{x}) = - \left( 4 \pi r \right) ^{-1}\) , we find: \({\mathbf {a}}(\mathbf{x}) = \nabla E(\mathbf{x})\), that is \({\mathbf {a}} \left( {\mathbf {x}} \right)\) is a dipole distribution, thus leading to Ladyzhenskaya’s solution:

$$\begin{aligned} p \left( {\mathbf {x}}, t \right) = {\mathbf {g}} \left( x, t \right) \cdot \widetilde{\mathbf{f}}_0, \end{aligned}$$
(38)

with

$$\begin{aligned} g_i(\mathbf{x},t) = - \frac{1}{4 \pi } \frac{\partial }{\partial x_i} \left( \frac{1}{r} \right) \delta \left( t \right) = \frac{x_i}{4 \pi r^3} \delta \left( t \right) \end{aligned}$$
(39)

3.1 Laplace domain solution

The gradient gauge can be applied on equal footing to the solution of the unsteady Stokes problem in the Laplace domain. Indicating with s the Laplace variable, and with a hat the Laplace transform of the corresponding quantity, Eq. (23) becomes

$$\begin{aligned} a^2 \widehat{\mathbf{v}}(\mathbf{x}) = \nabla ^2 \widehat{\mathbf{v}}(\mathbf{x}) - \frac{1}{\mu } \nabla {\widehat{p}}(\mathbf{x}) + \frac{\mathbf{f}_0}{\mu } \, \delta (\mathbf{x}-\mathbf{x}_0) \end{aligned}$$
(40)

where \(a=\sqrt{s/\nu }\), and the functional dependence on s as been omitted for notational simplicity. Equation (7) transforms into \({\widehat{v}}(\mathbf{x})={\widehat{v}}^{\, \prime }(\mathbf{x}) +\nabla {\widehat{\phi }}(\mathbf{x})\), where the freely evolving term \({\widehat{v}}^{\, \prime }(\mathbf{x})\) is a solution of the Helmholtz equation

$$\begin{aligned} \nabla ^2 {\widehat{v}}^{\, \prime }(\mathbf{x}) - a^2 \, {\widehat{v}}^{\, \prime }(\mathbf{x}) = -\frac{\mathbf{f}_0}{\mu } \, \delta (\mathbf{x}-\mathbf{x}_0) \end{aligned}$$
(41)

From the expression of the fundamental solution of the Helmholtz problem one obtains

$$\begin{aligned} {\widehat{v}}^{\, \prime }(\mathbf{x}) = \frac{e^{-a r}}{ 4 \pi \mu r} \, \mathbf{f}_0 \end{aligned}$$
(42)

where \(\mathbf{r}=\mathbf{x}-\mathbf{x}_0\). The gradient gauge field \({\widehat{\phi }}(\mathbf{x})\) is the solution of the equation

$$\begin{aligned} \nabla ^2 {\widehat{\phi }}(\mathbf{x}) = - \frac{1}{4 \pi \mu r} \frac{d }{d r} \left( \frac{e^{-ar}}{r} \right) \, (\mathbf{r} \cdot \mathbf{f}_0) \end{aligned}$$
(43)

that can be expressed as

$$\begin{aligned} {\widehat{\phi }}(\mathbf{x}) = \frac{1}{4 \pi \mu } u(r) \, (\mathbf{r} \cdot \mathbf{f}_0) \end{aligned}$$
(44)

Substituting the latter expression into Eq. (43), after simple but tedious calculations, the expression for u(r) follows

$$\begin{aligned} u(r)= -\frac{1}{a^2 \, r^3} + \frac{e^{-a r}}{a^2 \, r^3} + \frac{e^{-a r}}{a \, r^2} \end{aligned}$$
(45)

and the final expression for the Oseen tensor \({\widehat{G}}_{i,j}(\mathbf{x})\), \({\widehat{v}}_i(\mathbf{x})= \sum _{j} {\widehat{G}}_{i,j}(\mathbf{x}-\mathbf{x}_0) \, f_{0,j}\), is given by

$$\begin{aligned} {\widehat{G}}_{i,j}(\mathbf{x}) = \frac{1}{4 \pi \mu } \left[ \frac{e^{-ar}}{r} + u(r) \right] \, \delta _{i,j} + \frac{x_i x_j}{r} \frac{d u(r)}{d r} \end{aligned}$$
(46)

where, in Eq. (46) we have indicated \(r=|\mathbf{x}|\).

4 Discussion: infinite-propagation and the paradox of incompressibility

The results obtained in the previous section, coupled with the application of the gradient gauge, highlight the occurrence of a remarkable paradox in the hydrodynamic behavior of incompressible flows, which can be clearly appreciated from the regularity conditions at infinity satisfied by the solutions of the unsteady Stokes equation.

To clarify this issue, consider a scalar problem such as heat conduction described by the equation

$$\begin{aligned} \frac{\partial T(\mathbf{x},t)}{\partial t}= \alpha \, \nabla ^2 T(\mathbf{x},t) \end{aligned}$$

in \({{\mathbb {R}}}^3\), starting from an impulsive initial condition \(T(\mathbf{x},0)=T_0 \, \delta (\mathbf{x})\). Albeit, as for any parabolic models, the heat-conduction equation suffers of the problem of infinite propagation velocity [19, 20], its solutions satisfy nice regularity properties at infinity, namely that both \(T(\mathbf{x},t)\) and its gradient \(\nabla T(\mathbf{x},t)\) vanish, at infinity, faster than any power of \(|\mathbf{x}|\),

$$\begin{aligned} \lim _{|\mathbf{x}| \rightarrow \infty } |\mathbf{x}|^m \, T(\mathbf{x},t)= \lim _{|\mathbf{x}| \rightarrow \infty } |\mathbf{x}|^m \, | \nabla T(\mathbf{x},t) | =0 \end{aligned}$$
(47)

for any \(m=0,1,\dots\). Indeed, the solution of the above problem is given by \(T(\mathbf{x},t)=G_D(\mathbf{x},t) \, T_0\), where \(G_D(\mathbf{x},t)\) is expressed by Eq. (25), with \(\nu\) substituted by \(\alpha\). Consequently, while its parabolic nature poses intrinsic problems associated with its relativistic consistency, the regularity conditions (47) are sufficiently strong to support the application of the parabolic approximation for low-energy applications, i.e., in classical heat trasport problems of common engineering and physics interest.

Next, turn the attention to the corresponding hydrodynamic solution in the unsteady Stokes regime. In this case, the velocity field \(\mathbf{v}(\mathbf{x},t)\), solution of Eq. (23), does not possess the same regularity conditions at infinity observed in the scalar problem, but a much weaker condition, namely that \(\mathbf{v}(\mathbf{x},t)\) should vanish as \(|\mathbf{x}| \rightarrow \infty\). The explicit calculation of the Green function provides the large-scale property

$$\begin{aligned} |\mathbf{v}(\mathbf{x},t)| \sim \frac{1}{|\mathbf{x}|^3} \end{aligned}$$
(48)

for large \(|\mathbf{x}|\).

In the light of the manifold of hydrodynamic paradoxes involving Stokesian flows [18], such as those occurring in two-dimensional models for which a hydrodynamic solution consistent with the prescribed, and physically reasonable, boundary conditions does not exist, this result could be viewed as a mild inconsistency. At a more careful investigation, Eq. (48) determines an intrinsically paradoxical behavior of the incompressible time-dependent Stokes equation, in which the effect of unbounded propagation typical of parabolic models is amplified by the vectorial nature of the problem. This phenomenon can be clearly appreciated by adopting the decomposition (7), as the velocity field \(\mathbf{v}(\mathbf{x},t)\) is the superposition of the vector \(\mathbf{v}^\prime (\mathbf{x},t)\) that propagates as a scalar field fulfilling a parabolic equation (24), and of the gradient gauge \(\phi (\mathbf{x},t)\), that is the responsible for the “unphysical” scaling Eq. (48), associated with the solution of the elliptic equation (27). Indeed, the physical reason for this behavior is elementary: the gradient gauge, that is the basic ingredient in assessing incompressibility, does not propagate, it simply satisfies a Poisson equation, in which time enters as a mere parameter, and owing to the long-term scaling of its fundamental solution, determines the power-law scaling (48). Another view to the scaling relation (48) is presented in the next paragraph.

Therefore, while in scalar parabolic models, the paradox of infinite propagation velocity stems exclusively from the simplifications underlying the use of a Fickian constitutive equation for the fluxes proportional to the concentration gradients (determining the occurrence of a Laplacian operator in the balance equations, and thus the parabolic character of the evolution operator), in transport equations for vector fields (the velocity field) in the presence of dynamic constraints (such as incompressibility), the steady and timeless representation of these constraints, leads to an unphysical behavior in the large-distance limit.

This phenomenon can be further argumented by considering a simple example. Let us suppose to cure the unphysical unbounded propagation of viscous flows by generalizing the linearized Navier–Stokes equation by means of an operator \({{\mathcal {L}}}_\mathrm{fp}\) possessing finite propagation c, in the meaning that any initially compactly supported velocity field is mapped into a compactly supported velocity field at any later time \(t>0\). For example one can choose a Cattaneo-type equation [19,20,21] to describe this property

$$\begin{aligned} \rho \, \frac{\partial \mathbf{v}(\mathbf{x},t)}{\partial t} ={{\mathcal {L}}}_\mathrm{fp}[\mathbf{v}(\mathbf{x},t)] \end{aligned}$$
(49)

where

$$\begin{aligned} {{\mathcal {L}}}_\mathrm{fp}[\mathbf{v}(\mathbf{x},t)]&= \mu \, \nabla ^2 \mathbf{v}(\mathbf{x},t) -\rho \, \tau _c \, \frac{\partial ^2 \mathbf{v}(\mathbf{x},t)}{\partial t^2} -\nabla p(\mathbf{x},t) \nonumber \\&\quad+ \mathbf{f}_0 \, \delta (\mathbf{x}-\mathbf{x}_0) \, \delta (t) \end{aligned}$$
(50)

and \(\tau _c\) is some constant characteristic time. It is important to observe that, albeit the Cattaneo model possesses some intrinsic inconsistencies for spatial dimensions greater than one (lack of positivity) [22, 23], it can fit the purposes of the present analysis in which solely the condition of finite propagation velocity is required.

Therefore,   the hydrodynamic problem considered consists in the Eq. (50), equipped with the incompressibility condition, and with vanishing initial conditions. To Eq. (50) the gauge decomposition (7) can be applied, and the \(\mathbf{v}^\prime\)-component satisfies the scalar equation

$$\begin{aligned} \rho \, \left[ \frac{\partial \mathbf{v}^\prime (\mathbf{x},t)}{\partial t} \right.+ & {} \left. \tau _c \, \frac{\partial ^2 \mathbf{v}^\prime (\mathbf{x},t)}{\partial t^2} \right] = \mu \, \nabla ^2 \mathbf{v}^\prime (\mathbf{x},t) \nonumber \\+ & {} \mathbf{f}_0 \, \delta (\mathbf{x}-\mathbf{x}_0) \, \delta (t) \end{aligned}$$
(51)

The corresponding solution is thus given by

$$\begin{aligned} \mathbf{v}^\prime (\mathbf{x},t) = G_\mathrm{Cat}(\mathbf{x}-\mathbf{x}_0,t) \, \widetilde{\mathbf{f}}_0 \end{aligned}$$
(52)

where \(G_\mathrm{Cat}(\mathbf{x},t)\) is the scalar Green function for the three-dimensional Cattaneo equation, that can be found e.g. in [24], which is compactly supported and vanishing, in the present case, outside the ball \(B(R(t),\mathbf{x}_0)\) of radius R(t) centered at \(\mathbf{x}_0\), where

$$\begin{aligned} R(t)= c \, t \, , \qquad c= \sqrt{\frac{\nu }{\tau _c}} \end{aligned}$$
(53)

Next, consider the gradient gauge \(\phi (\mathbf{x},t)\), that is the solution of the equation deriving from incompressibility

$$\begin{aligned} \nabla ^2 \phi (\mathbf{x},t)= - \nabla _x G_\mathrm{Cat}(\mathbf{x}-\mathbf{x}_0,t) \cdot \widetilde{\mathbf{f}}_0 \end{aligned}$$
(54)

and thus

$$\begin{aligned} \phi (\mathbf{x},t)= \frac{1}{4 \, \pi } \int _{B(R(t),\mathbf{x}_0)} \frac{\nabla _{\xi } G_\mathrm{Cat}(\varvec{\xi }-\mathbf{x}_0,t) \cdot \widetilde{\mathbf{f}}_0}{|\mathbf{x}-\varvec{\xi }|} \, d \varvec{\xi } \end{aligned}$$
(55)

Independently of the fact that the viscous propagation described by Eq. (51) possesses finite propagation velocity, the gauge field scales for large \(|\mathbf{x}-\mathbf{x}_0|\) as \(\phi (\mathbf{x}) \sim 1/|\mathbf{x}-\mathbf{x}_0|^2\), and since \(\mathbf{v}(\mathbf{x},t)=\mathbf{v}^\prime (\mathbf{x},t)+ \nabla \phi (\mathbf{x},t)\), Eq. (48) follows for large \(|\mathbf{x}|\) and for any \(t>0\). This result is a direct consequence of the gradient gauge, forcing the scalar field \(\phi (\mathbf{x},t)\) to satisfy an elliptic problem. It indicates that the power-law tail characterizing the Green function of linear unsteady hydrodynamic problems is related to the incompressibility assumption that generates non-local effects and unbounded propagation even in the presence of hyperbolic operators for describing the internal viscous stresses. In other words, any physically consistent model providing the bounded propagation of the viscous stresses necessarily implies the removal of the unphysical condition of incompressibility, which can be viewed as an equilibrium approximation, i.e. valid for liquids at rest. Said in different words, even for liquids, the large-distance propagation of any disturbance no matter how small it is, should necessarily imply compressibility effects, which does not mean that the fluid should behave as inviscid, but solely that the assumption of a solenoidal velocity field is not compatible with the finite propagation of the internal disturbances. This observation is relevant in the hydrodynamic theory of Brownian motion, and in the analysis of transient properties of colloidal suspensions.

The hydrodynamic theory of Brownian motion [25, 26] extended the original approach due to Einstein and Langevin [27, 28] based on the instantaneous Stokes equation, by considering fluid-particle interactions in the time-dependent incompressible Stokes regime, thus including the inertial effects in the fluid. As a result of this improvement, in has been shown that the velocity autocorrelation of a Brownian particle deviates from the exponential relaxation predicted by the Einstein-Langevin model, displaying long-term power law decay with time characterized, for Newtonian fluids, by the exponent 3/2, deriving from the presence of the Basset force. Moreover, as regards the mean square deviation of velocity fluctuations, the incompressible hydrodynamic theory predicts, at constant temperature T, the value \(k_B \, T/(m+m_a)\), where m is the particle mass, \(m_a\) the hydrodynamic added mass, and \(k_B\) the Boltzmann constant, providing a qualitative different prediction from the classical result \(k_B \, T/m\) of classical statistical physics (law of energy equipartition). In order to resolve this contradiction, Zwanzing and Bixon [29, 30] and Chow and Hermans [31] showed that if incompressibility is removed, by assuming an equation of state for pressure proportial to the density (through the square of the velocity of sound), the classical statistical physical predictions are recovered. Recent experimental works on the fine structure of Brownian fluctuations, resolving time scales below the characteristic dissipation time (order of \(10^{-7}\) s for micrometric particles in water at room temperature) have shown not only the quantitative agreement with the hydrodynamic theory of Brownian motion as regards the temporal behavior of the velocity autocorrelation function, but also the violation of energy equipartition predicted by the added-mass effect [32,33,34]. These experiments and their physical relevance suggest for an improvement of the present description of low-velocity hydrodynamics at microscale, in which not only the incompressibility assumption is removed, but also such that it should be able to predict fluid-particle interactions in any experimental conditions involving pressure-driven flows in microchannels, still satisfying the physical requirement of a correct description of the acoustic propagation. The existing weak-compressibility models, essentially based on the perturbative expansion of the acoustic contributions around an incompressible leading-order term [35] do not fulfil this program. A discussion on this issue is posponed to paragraph 4.2.

4.1 Helmholtz decomposition of the vectorial Dirac’s \(\delta\)

There is a simple and striking proof that the Stokesian paradox associated with incompressibility is a consequence of the vectorial nature of the velocity field, owing to the Helmholtz decomposition in \({{\mathbb {R}}}^3\). Consider again the unsteady Stokes equation (24) in the presence of an impulsive forcing, setting \(\mathbf{x}_0={{\mathbf {0}}}\) for notational simplicity. To this forcing field, Helmholtz decomposition can be applied, i.e.,

$$\begin{aligned} \mathbf{f}_0 \, \delta (\mathbf{x}) = - \nabla \psi _\delta (\mathbf{x}) + \nabla \times \mathbf{A}_\delta (\mathbf{x}) \end{aligned}$$
(56)

where \(\psi _\delta (\mathbf{x})\) and \(\mathbf{A}_\delta (\mathbf{x})\) are the scalar and vector potentials of a vectorial Dirac’s \(\delta\) of intensity \(|\mathbf{f}_0|\) and oriented in the direction of \(\mathbf{f}_0\),

$$\begin{aligned} \psi _\delta (\mathbf{x})= & {} \frac{1}{4 \pi } \int _{{{\mathbb {R}}}^3} \frac{\nabla ^\prime \cdot \left( \mathbf{f}_0 \delta (\mathbf{x}^\prime ) \right) }{|\mathbf{x}-\mathbf{x}^\prime |} \, d \mathbf{x}^\prime \nonumber \\ \mathbf{A}_\delta (\mathbf{x})= & {} \frac{1}{4 \pi } \int _{{{\mathbb {R}}}^3} \frac{\nabla ^\prime \times \left( \mathbf{f}_0 \delta (\mathbf{x}^\prime ) \right) }{|\mathbf{x}-\mathbf{x}^\prime |} \, d \mathbf{x}^\prime \end{aligned}$$
(57)

where \(\nabla ^\prime\) indicates the nabla operator acting on the \(\mathbf{x}^\prime\)-coordinates. To begin with, consider \(\psi _\delta (\mathbf{x})\). Differentiating by parts the integrand at the r.h.s. of the Helmholtz representation of \(\psi _\delta (\mathbf{x})\),

$$\begin{aligned} \psi _\delta (\mathbf{x})= & {} \frac{1}{4 \pi } \int _{{{\mathbb {R}}}^3} \nabla ^\prime \cdot \left( \frac{ \mathbf{f}_0 \delta (\mathbf{x}^\prime )}{|\mathbf{x}-\mathbf{x}^\prime |} \right) \, d \mathbf{x}^\prime \nonumber \\- & {} \frac{1}{4 \pi } \int _{{{\mathbb {R}}}^3} \frac{\delta (\mathbf{x}^\prime )}{|\mathbf{x}-\mathbf{x}^\prime |^2} \, \mathbf{f}_0 \cdot \nabla (|\mathbf{x}-\mathbf{x}^\prime |) \, d \mathbf{x}^\prime \end{aligned}$$
(58)

in which the property \(\nabla h(\xi )=-\nabla ^\prime h(\xi )\) for any function \(h(\xi )\) of argument \(\xi =|\mathbf{x}-\mathbf{x}^\prime |\) has been applied. The first integral at the r.h.s. of Eq. (58) can be transformed by the Gauss divergence theorem into a surface integral on the spherical surface at infinity, and is identically vanishing, owing to the localization properties of \(\delta (\mathbf{x}^\prime )\). Thus Eq. (58) simplifies as

$$\begin{aligned} \psi _\delta (\mathbf{x}) =- \frac{1}{4 \pi \, |\mathbf{x}|^2} \mathbf{f}_0 \cdot \nabla (|\mathbf{x}|)=- \frac{\mathbf{x} \cdot \mathbf{f}_0}{4 \pi |\mathbf{x} |^3} \end{aligned}$$
(59)

Analogously for the vector potential, reworking by parts the integrand, one obtains

$$\begin{aligned} \mathbf{A}_\delta (\mathbf{x})= & {} \frac{1}{4 \pi } \int _{{{\mathbb {R}}}^3} \nabla ^\prime \times \left( \frac{ \mathbf{f}_0 \delta (\mathbf{x}^\prime )}{|\mathbf{x}-\mathbf{x}^\prime |} \right) \, d \mathbf{x}^\prime \nonumber \\- & {} \frac{1}{4 \pi } \int _{{{\mathbb {R}}}^3} \frac{\delta (\mathbf{x}^\prime )}{|\mathbf{x}-\mathbf{x}^\prime |^2} \, \nabla (|\mathbf{x}-\mathbf{x}^\prime |) \times \mathbf{f}_0\, d \mathbf{x}^\prime \end{aligned}$$
(60)

The first integral at the r.h.s of Eq. (60) is identically vanishing, as it can be transformed using the curl-version of the Gauss theorem into a surface integral on the spherical surface at infinity. Consequently, the second term at the r.h.s. of Eq. (60) yields

$$\begin{aligned} \mathbf{A}_\delta (\mathbf{x}) = - \frac{1}{4 \pi \, |\mathbf{x}|^2} \nabla (|\mathbf{x}|) \times \mathbf{f}_0 =- \frac{\mathbf{x} \times \mathbf{f}_0}{4 \pi |\mathbf{x} |^3} \end{aligned}$$
(61)

The physical implications of the decomposition (56) and of the expressions for \(\psi _\delta (\mathbf{x})\) and \(\mathbf{A}_\delta (\mathbf{x})\) given by Eqs. (59) and (61) are evident, by considering that the functional space \(L^2_{\mathrm{vec}}( {{\mathbb {R}}}^3)\) of vectorial square summable functions in \({{\mathbb {R}}}^3\) is the direct sum \(L^2_{\mathrm{vec}}( {{\mathbb {R}}}^3) = L^2_\mathrm{sol}( {{\mathbb {R}}}^3) \oplus L^2_\mathrm{irr}( {{\mathbb {R}}}^3)\), where

$$\begin{aligned} L^2_\mathrm{sol}( {{\mathbb {R}}}^3)= & {} \left\{ \, \mathbf{f}(\mathbf{x}) \; \left| \, \int _{{{\mathbb {R}}}^3} |\mathbf{f}(\mathbf{x}) |^2 \, d \mathbf{x} \; , \nabla \cdot \mathbf{f}(\mathbf{x})=0 \right. \right\} \nonumber \\ L^2_\mathrm{irr}( {{\mathbb {R}}}^3)= & {} \left\{ \, \mathbf{f}(\mathbf{x}) \; \left| \, \int _{{{\mathbb {R}}}^3} |\mathbf{f}(\mathbf{x}) |^2 \, d \mathbf{x} \; , \nabla \times \mathbf{f}(\mathbf{x})=0 \, , \right. \right. \nonumber \\&\,&\nabla \cdot \mathbf{f}(\mathbf{x}) \ne 0 \Big \} \end{aligned}$$
(62)

and that \(L^2_\mathrm{sol}( {{\mathbb {R}}}^3) \perp L^2_\mathrm{irr} ( {{\mathbb {R}}}^3)\), in the meaning that if \(\mathbf{v_s}(\mathbf{x}) \in L^2_\mathrm{sol}( {{\mathbb {R}}}^3)\) and \(\mathbf{v_i}(\mathbf{x}) \in L^2_\mathrm{irr}( {{\mathbb {R}}}^3)\), then \(\int _{{{\mathbb {R}}}^3} \mathbf{v}_s(\mathbf{x}) \cdot \mathbf{v}_i(\mathbf{x}) \, d \mathbf{x} =0\) [4].

The decomposition (56), substituted into Eq. (24), provides (with \(\mathbf{x}_0={{\mathbf {0}}}\)),

$$\begin{aligned} \rho \, \frac{\partial \mathbf{v}(\mathbf{x},t)}{\partial t}& = \mu \, \nabla ^2 \mathbf{v}(\mathbf{x},t) - \nabla p(\mathbf{x},t) \nonumber \\&\quad+ \left[ - \nabla \psi _\delta (\mathbf{x}) + \nabla \times \mathbf{A}_\delta (\mathbf{x}) \right] \, \delta (t) \end{aligned}$$
(63)

that, owing to the incompressibility of \(\mathbf{v}(\mathbf{x})\), and to the orthogonality of \(L^2_\mathrm{sol}( {{\mathbb {R}}}^3)\) and \(L^2_\mathrm{irr}( {{\mathbb {R}}}^3)\), decouples into two equations, one for each subspace of \(L^2_{\mathrm{vec}}( {{\mathbb {R}}}^3)\): (i) in \(L^2_\mathrm{irr}( {{\mathbb {R}}}^3)\), it provides direcly the equation for pressure (modulo an additive constant),

$$\begin{aligned} p(\mathbf{x},t)=- \psi _\delta (\mathbf{x}) \delta (t) \end{aligned}$$
(64)

coinciding with Eq. (39), (ii) in \(L^2_\mathrm{sol}( {{\mathbb {R}}}^3)\), it yields the evolution equation of the divergence-free velocity field

$$\begin{aligned} \rho \, \frac{\partial \mathbf{v}(\mathbf{x},t)}{\partial t} = -\mu \, \nabla \times \nabla \times \mathbf{v}(\mathbf{x},t) + \nabla \times \mathbf{A}_\delta (\mathbf{x}) \delta (t) \end{aligned}$$
(65)

It follows from Eqs. (56) and (64) that in the unsteady Stokes problem with impulsive forcing, the pressure should be also impulsive in time, corresponding to an instantaneous acoustic propagation (ipso facto a kind of “teleportation”). If pressure, as it is, is a physical variable related to the propagation of material fluid elements, the impulsive response of \(p(\mathbf{x},t)\), that follows from the geometry of \(L^2_{\mathrm{vec}}({{\mathbb {R}}}^3)\) is the most striking manifestation of the teleportation paradox implied by incompressibility. Equation (65) indicates that an impulsive perturbation, once projected onto \(L^2_\mathrm{sol}( {{\mathbb {R}}}^3)\) provides an effective forcing decaying with \(|\mathbf{x}|\) as \(1/|\mathbf{x}|^2\), and this is the geometric reason for the anomalous long-term scaling expressed by Eq. (48).

To conclude, consider the Helmholtz decomposition of the vector field \(\mathbf{v}^\prime (\mathbf{x},t)\) solution of Eq. (24): \(\mathbf{v}^\prime (\mathbf{x},t)=\mathbf{v}^\prime _s(\mathbf{x},t)+\mathbf{v}^\prime _i(\mathbf{x},t)\), where \(\mathbf{v}^\prime _s(\mathbf{x},t) \in L^2_\mathrm{sol}( {{\mathbb {R}}}^3)\), and \(\mathbf{v}^\prime _i(\mathbf{x},t) \in L^2_\mathrm{irr}( {{\mathbb {R}}}^3)\). Because of the mutual orthogonality of these subspaces, the problem is split into two separate equations

$$\begin{aligned} \rho \, \frac{\partial \mathbf{v}_s^\prime (\mathbf{x},t)}{\partial t}= - \mu \, \nabla \times \nabla \times \mathbf{v}^\prime _s(\mathbf{x},t) + \nabla \times \mathbf{A}_\delta (\mathbf{x}) \, \delta (t) \end{aligned}$$
(66)

and

$$\begin{aligned} \rho \, \frac{\partial \mathbf{v}_i^\prime (\mathbf{x},t)}{\partial t}= \mu \, \nabla ^2 \mathbf{v}_i^\prime (\mathbf{x},t) - \nabla \psi _\delta (\mathbf{x}) \, \delta (t) \end{aligned}$$
(67)

Since \(\mathbf{v}_s^\prime (\mathbf{x},t) = \nabla \times \mathbf{A}_v(\mathbf{x},t)\), \(\mathbf{v}_i^\prime (\mathbf{x},t)= \nabla \phi _v(\mathbf{x},t)\), the vector and scalar velocity potentials \(\mathbf{A}_v(\mathbf{x},t)\) and \(\phi _v(\mathbf{x},t)\) satisfy the equations

$$\begin{aligned} \rho \, \frac{\partial \mathbf{A}_v(\mathbf{x},t)}{\partial t}&= \mu \, \nabla ^2 \mathbf{A}_v(\mathbf{x},t) + \mathbf{A}_\delta (\mathbf{x},t) \, \delta (t) \nonumber \\ \rho \, \frac{\partial \phi _v(\mathbf{x},t)}{\partial t}&= \mu \, \nabla ^2 \phi _v(\mathbf{x},t) - \psi _\delta (\mathbf{x}) \, \delta (t) \end{aligned}$$
(68)

From the functional form of the forcing terms \(\mathbf{A}_\delta (\mathbf{x})\) and \(-\psi _\delta (\mathbf{x})\), see Eqs. (59) and (61), corresponding to the projection in the two Helmholtz subspaces of an impulsive vectorial perturbation, both \(\mathbf{v}^\prime _s(\mathbf{x},t)\) and \(\mathbf{v}^\prime _i(\mathbf{x},t)\) are characterized, at any time \(t>0\), by long-distance power-law tails Eq. (48), that mutually cancel out once their superposition \(\mathbf{v}^\prime (\mathbf{x},t)\) is considered (providing in the large-distance limit the typical scaling of a heat kernel, i.e. \(\mathbf{v}^\prime (\mathbf{x},t) \sim e^{-|\mathbf{x}-\mathbf{x}_0|^2/4 \, \nu \, t}\)). It can be concluded that the Helmholtz decomposition itself generates at any time \(t>0\) large-distance effects as regards the projection of \(\mathbf{v}^\prime (\mathbf{x},t)\) onto the subsets \(L^2_\mathrm{sol}( {{\mathbb {R}}}^3)\) and \(L^2_\mathrm{irr}( {{\mathbb {R}}}^3)\), and these effects vanish once \(\mathbf{v}^\prime (\mathbf{x},t)\) is considered in its unconstrained evolution.

4.2 Gradient-gauge decomposition and the role of pressure

The gradient-gauge formulation of incompressible Stokesian dynamics is essentially aimed at showing the inescapable paradoxes in the instantaneous propagation of hydrodynamic fields associated with the assumption of incompressibility. But the origin of the gradient-gauge formulation stems from the ancillary role of the pressure in incompressible hydrodynamics, that albeit in a slightly different mathematical setting, plays the role of a gradient-gauge variable, just to enforce incompressibility. The existing pertubative expansions of weak-compressible theory [35] do not resolve the above mentioned paradoxes in the most general setting of an isothermal pressure-driven flow interacting with a micrometric particle, as shown by the recent literature on Brownian motion fluctuations. A more accurate description of fluid-particle interactions at the microscale is thus required in the analysis of hydrodynamic processes below the dissipation timescale (order of \(10^{-7}\) s, for a micrometric particle in water at room temperature).

The analysis of the Helmholtz decomposition of a vectorial impulsive perturbation addressed in the previous paragraph, and of its consequences as regards the large-distance scaling of the two orthogonal contributions lying in the two subspaces \(L^2_\mathrm{sol}({{\mathbb {R}}}^3)\) and \(L^2_\mathrm{irr}({{\mathbb {R}}}^3)\), paves the way towards the improvement of the actual hydrodynamic formulation of low-velocity flows able to account for the phenomenologies occurring at microscale (Brownian motion and motion of particles in microchannels).

The first prerequisite is that the unphysical incompressibility assumption should be removed, permitting a correct description of acoustic perturbations below the dissipation timescale. From the analysis developed at the end of paragraph 4.1, such formulation should involve the free propagation of the velocity field, not subjected to any constraints associated with the Helmholtz decomposition. A way for achieving this program, without altering the structure of the Stokes or Navier–Stokes equations (that provide an excellent description of macroscopic flows) is: (i) to restore to pressure its full nature of physical variable accounting for the internal compression stresses in a liquid far from mechanical equilibrium, and (ii) to derive for it an evolution equation that, in the limit case of steady flows should recover incompressibility as a limit case. This approach is outlined in [36] and does not lead neither to conceptual or structural modifications of the Navier–Stokes equations, but rather to their completion via the introduction of a physically-derived evolution dynamics for the pressure variable.

5 Hydrodynamic Green functions in the presence of boundaries

Apart from the analysis of the paradoxes associated with incompressible unsteady Stokes flows, the gradient gauge can be used to obtain a formal solution (in the form of an integral equation) for the hydrodynamic Green functions in the presence of solid boundaries. Henceforth, no-slip boundary conditions are assumed. This section addresses the mathematical structure of the integral equations for the hydrodynamic Green function in Stokes and time-dependent Stokes regimes, emerging from the gradient gauge approach, for flows defined on \(\Omega \subseteq {{\mathbb {R}}}^n\), where \(n=2,3\), in the case \(\Omega\) possesses a solid boundary \(\partial \Omega\) at which the velocity vanishes. In virtue of the discussion of the previous section, in order to avoid the classical Stokesian paradoxes [18], the two-dimensional case, \(n=2\), is limited to bounded domains \(\Omega\).

To begin with, consider the Stokes problem for the hydrodynamic Green functions, defined by Eqs. (5–6) for \(\mathbf{x}, \mathbf{x}_0 \in \Omega\), and such that for any \(\varvec{\xi } \in \partial \Omega\),

$$\begin{aligned} \mathbf{v}(\varvec{\xi }) = 0 \, , \qquad \varvec{\xi } \in \partial \Omega \end{aligned}$$
(69)

In the presence of solid boundaries, the decomposition (7) valid for the free-space propagation should be generalized as

$$\begin{aligned} \mathbf{v}(\mathbf{x})=\mathbf{v}^\prime (\mathbf{x}) + \mathbf{v}^{\prime \prime }(\mathbf{x}) + \nabla \phi (\mathbf{x}) \end{aligned}$$
(70)

where \(\mathbf{v}^\prime (\mathbf{x})\) is the solution of the elliptic problem Eq. (8), defined in \(\Omega\) with the boundary condition at \(\partial \Omega\),

$$\begin{aligned} \mathbf{v}^\prime (\varvec{\xi }) =0 \, \qquad \varvec{\xi } \in \partial \Omega \end{aligned}$$
(71)

and \(\mathbf{v}^{\prime \prime }(\mathbf{x})\) is the solution of the vector-valued Laplace equation in \(\Omega\)

$$\begin{aligned} \nabla ^2 \mathbf{v}^{\prime \prime }(\mathbf{x})=0 \end{aligned}$$
(72)

equipped with generic Dirichlet boundary conditions,

$$\begin{aligned} \mathbf{v}^{\prime \prime }(\varvec{\xi })= \mathbf{v}_b({\varvec{\xi }}) \, , \qquad \varvec{\xi } \in \partial \Omega \end{aligned}$$
(73)

The boundary function \(\mathbf{v}_b(\varvec{\xi })\), is the extra degree of freedom introduced in order to match the no-slip boundary conditions at \(\partial \Omega\), and it will be determined subsequently. From Eqs. (8) and (72) it follows that the gradient gauge \(\phi (\mathbf{x})\) satisfies in \(\Omega\) the Poisson equation

$$\begin{aligned} \nabla ^2 \phi (\mathbf{x}) = - \nabla \cdot \mathbf{v}^\prime (\mathbf{x}) - \nabla \cdot \mathbf{v}^{\prime \prime }(\mathbf{x}) \end{aligned}$$
(74)

For this equation, assume at \(\partial \Omega\) homogeneous Dirichlet conditions

$$\begin{aligned} \phi (\varvec{\xi })=0 \, , \qquad \ \varvec{\xi } \in \partial \Omega \end{aligned}$$
(75)

Also in this case the pressure \(p(\mathbf{x})\) is defined by Eq. (10). Consider the boundary condition at \(\partial \Omega\) induced by the decomposition (70). Because of Eqs. (71), (73), it reduced to

$$\begin{aligned} \mathbf{v}^{\prime \prime }(\varvec{\xi }) + \nabla \phi (\mathbf{x}) |_{\mathbf{x}=\varvec{\xi }} =0 \, , \qquad \varvec{\xi } \in \partial \Omega \end{aligned}$$
(76)

Let \(\mathbf{t}(\varvec{\xi })\) be any tangent unit vector on the boundary manifold \(\partial \Omega\), and \(\mathbf{n}_e(\varvec{\xi })\) the outer normal unit vector at \(\varvec{\xi } \in \partial \Omega\). Because of Eq. (75), \(\mathbf{t}(\varvec{\xi }) \cdot \nabla \phi (\mathbf{x}) |_{\mathbf{x}=\varvec{\xi }}=0\) identically, so that Eq. (76) reduces to

$$\begin{aligned} \mathbf{v}_b(\varvec{\xi }) \cdot \mathbf{t}(\varvec{\xi }) = 0 \end{aligned}$$
(77)
$$\begin{aligned} \mathbf{v}_b(\varvec{\xi }) \cdot \mathbf{n}_e(\varvec{\xi }) + \left. \frac{\partial \phi (\mathbf{x})}{\partial n_e} \right| _{\mathbf{x} = \varvec{\xi }} = 0 \end{aligned}$$
(78)

for any \(\varvec{\xi } \in \partial \Omega\). Condition (77) indicates that \(\mathbf{v}_b(\varvec{\xi })\) is purely normal to \(\partial \Omega\). Consequently, it is sufficient to introduce the scalar function \(\psi (\varvec{\xi })\) at \(\partial \Omega\), such that

$$\begin{aligned} \mathbf{v}_b(\varvec{\xi })= \psi (\varvec{\xi }) \, \mathbf{n}_e(\varvec{ \xi }) \end{aligned}$$
(79)

This auxiliary boundary function \(\psi (\varvec{\xi })\) is determined by the remaining boundary condition (78), dictating

$$\begin{aligned} \psi (\varvec{\xi }) \left. \frac{\partial \phi (\mathbf{x})}{\partial n_e} \right| _{\mathbf{x} = \varvec{\xi }} = 0 \, , \qquad \varvec{\xi } \in \partial \Omega \end{aligned}$$
(80)

The mathematical structure of the problem is set, and its solution involves solely the fundamental solutions of the scalar Poisson and Laplace equations in \(\Omega\). More precisely, let \(G(\mathbf{x},\mathbf{x}_0)\) the Green function for the scalar elliptic problem

$$\begin{aligned} \nabla ^2 G(\mathbf{x},\mathbf{x}_0) = - \delta (\mathbf{x}-\mathbf{x}_0) \end{aligned}$$
(81)

equipped with the homogeneous Dirichlet condition

$$\begin{aligned} G(\varvec{\xi },\mathbf{x_0}) =0 \, , \qquad \varvec{\xi } \in \partial \Omega \end{aligned}$$
(82)

and introduce further the fundamental solution \(E(\mathbf{x},\varvec{\xi })\) for the scalar Laplace equation in \(\Omega\)

$$\begin{aligned} \nabla ^2 E(\mathbf{x},\varvec{\xi })&= 0 \, , \qquad \mathbf{x} \in \Omega \, , \;\;\; \varvec{\xi } \in \partial \Omega \nonumber \\ E(\varvec{\eta },\varvec{\xi })&= \delta (\varvec{\eta }-\varvec{\xi }) \, , \qquad \varvec{\eta } \in \partial \Omega \end{aligned}$$
(83)

Next, consider the determination of \(\mathbf{v}^\prime (\mathbf{x})\), \(\mathbf{v}^{\prime ,\prime }(\mathbf{x})\) and \(\phi (\mathbf{x})\).

From Eq. (8) and (81–82) it follows that \(\mathbf{v}^\prime (\mathbf{x})\) is given by

$$\begin{aligned} \mathbf{v}^\prime (\mathbf{x})= \frac{G(\mathbf{x},\mathbf{x}_0) \, \mathbf{f}_0}{\mu } \end{aligned}$$
(84)

As regards \(\mathbf{v}^{\prime \prime }(\mathbf{x})\), Eqs. (79) and (83) provide

$$\begin{aligned} \mathbf{v}^{\prime \prime }(\mathbf{x}) = \int _{\partial \Omega } E(\mathbf{x},\varvec{\xi }) \, \psi (\varvec{\xi }) \, \mathbf{n}_e (\varvec{\xi }) \, d S(\varvec{\xi }) \end{aligned}$$
(85)

where \(d S(\varvec{\xi })\) indicates the surface element on \(\partial \Omega\). Finally, the formal solution of the gradient gauge Eq. (74) is given by

$$\begin{aligned} \phi (\mathbf{x})= \int _{\Omega } G(\mathbf{x},\mathbf{y}) \, \nabla _y \cdot \mathbf{v}^\prime (\mathbf{y}) \, d \mathbf{y} +\int _{\Omega } G(\mathbf{x},\mathbf{y}) \, \nabla _y \cdot \mathbf{v}^{\prime \prime }(\mathbf{y}) \, d \mathbf{y} \end{aligned}$$
(86)

It remains to determine \(\psi (\varvec{\xi })\). Enforcing the representations (84–86) within the boundary condition (80), one finally obtains for \(\psi (\varvec{\xi })\) the Fredholm integral equation of the second kind

$$\begin{aligned} \psi (\varvec{\xi }) + m(\varvec{\xi }) + \int _{\partial \Omega } K(\varvec{\xi },\varvec{\eta }) \, \psi (\varvec{\eta }) \, d S(\varvec{\eta }) =0 \end{aligned}$$
(87)

where

$$\begin{aligned} m(\varvec{\xi })&= \frac{1}{\mu } \int _{\Omega } \mathbf{n}_e(\varvec{\xi }) \cdot \nabla _x G(\mathbf{x},\mathbf{y} )|_{\mathbf{x}=\varvec{\xi }} \, \nabla _y G(\mathbf{y}, \mathbf{x}_0) \cdot \mathbf{f}_0 \, d \mathbf{y} \end{aligned}$$
(88)
$$\begin{aligned} K(\varvec{\xi },\varvec{\eta })&= \int _{\Omega } \mathbf{n}_e(\varvec{\xi }) \cdot \nabla _x G(\mathbf{x},\mathbf{y} )|_{\mathbf{x}=\varvec{\xi }} \, \nabla _y E(\mathbf{y}, \varvec{\eta }) \cdot \mathbf{n}_e(\varvec{\eta }) \, d \mathbf{y} \end{aligned}$$
(89)

The same approach can be transferred to the analysis of the time-dependent Stokes problem Eq. (23) on \(\Omega\), where in the present case the boundary condition on \(\partial \Omega\) reads as

$$\begin{aligned} \mathbf{v}(\varvec{\xi },t) =0 \, , \qquad \varvec{\xi } \in \partial \Omega \end{aligned}$$
(90)

and \(\mathbf{v}(\mathbf{x},0)=0\). Enforcing the same decomposition, introduced for the steady Stokes problem,

$$\begin{aligned} \mathbf{v}(\mathbf{x},t)=\mathbf{v}^\prime (\mathbf{x},t) + \mathbf{v}^{\prime \prime }(\mathbf{x},t) + \nabla \phi (\mathbf{x},t) \end{aligned}$$
(91)

the component \(\mathbf{v}^\prime (\mathbf{x},t)\) is the solution of Eq. (24) with \(\mathbf{v}^\prime (\varvec{\xi },t)=0\) for \(\varvec{\xi } \in \partial \Omega\) and \(\mathbf{v}^\prime (\mathbf{x},0)=0\), while \(\mathbf{v}^{\prime \prime }(\mathbf{x},t)\) is the solution of the homogeneous vector-valued diffusion equation

$$\begin{aligned} \frac{\partial \mathbf{v}^{\prime \prime }(\mathbf{x},t)}{\partial t} = \nu \, \nabla ^2 \mathbf{v}^{\prime \prime }(\mathbf{x},t) \end{aligned}$$
(92)

equipped with vanishing initial conditions \(\mathbf{v}^{\prime \prime }(\mathbf{x},0)=0\) and with the Dirichlet boundary condition

$$\begin{aligned} \mathbf{v}^{\prime \prime }(\varvec{\xi },t)= \psi (\varvec{\xi },t) \, \mathbf{n}_e(\varvec{\xi }) \, , \qquad \varvec{\xi } \in \partial \Omega \end{aligned}$$
(93)

The gradient gauge \(\phi (\mathbf{x},t)\) satisfies Eq. (74), where all the fields entering this equation explicitly depend on time t, equipped with the homogeneous Dirichet condition at \(\partial \Omega\).

As for steady Stokes flow, the scalar boundary function \(\psi (\varvec{\xi },t)\) satisfies the boundary condition (80) for \(t>0\), replacing \(\psi (\varvec{\xi })\) and \(\phi (\mathbf{x})\) with \(\psi (\varvec{\xi },t)\) and \(\phi (\mathbf{x},t)\), respectively.

Let us indicate with \(G_D(\mathbf{x},\mathbf{x}_0,t)\) the Green function for the parabolic scalar problem

$$\begin{aligned} \frac{\partial G_D(\mathbf{x},\mathbf{x}_0,t)}{\partial t} = \nu \, \nabla _x^2 G_D(\mathbf{x},\mathbf{x}_0,t)+ \delta (\mathbf{x}-\mathbf{x}_0) \, \delta (t) \end{aligned}$$
(94)

with \(G_D(\mathbf{x},\mathbf{x}_0,0)=0\) and \(G_D(\varvec{\xi },\mathbf{x}_0,t)=0\) for \(\varvec{\xi } \in \partial \Omega\), and with \(E_t(\mathbf{x},\varvec{\xi },t)\) the fundamental solution of the parabolic scalar problem

$$\begin{aligned} \frac{\partial E_t(\mathbf{x},\varvec{\xi },t)}{\partial t}&= \nu \, \nabla _x^2 E_t(\mathbf{x},\varvec{\xi },t) \qquad \mathbf{x} \in \Omega \, , \;\;\; \varvec{\xi } \in \partial \Omega \nonumber \\ E_t(\varvec{\eta },\varvec{\xi },t)&= \delta (\varvec{\eta }- \varvec{\xi }) \delta (t) \, , \qquad \varvec{\eta } \in \partial \Omega \end{aligned}$$
(95)

and \(E_t(\mathbf{x},\varvec{\xi },0)=0\). Using these functions, the velocity field \(\mathbf{v}^\prime (\mathbf{x},t)\) can be expressed as

$$\begin{aligned} \mathbf{v}^\prime (\mathbf{x},t)= G_D(\mathbf{x},\mathbf{x}_0,t) \, \widetilde{\mathbf{f}}_0 \end{aligned}$$
(96)

with \(\widetilde{\mathbf{f}}_0=\mathbf{f}_0/\rho\), and

$$\begin{aligned} \mathbf{v}^{\prime \prime }(\mathbf{x},t)= \int _0^t d \tau \int _{\partial \Omega } E_t(\mathbf{x},\varvec{\xi },t-\tau ) \, \psi (\varvec{\xi },\tau ) \, \mathbf{n}_e(\varvec{\xi }) \, d S(\varvec{\xi }) \end{aligned}$$
(97)

Consequently,

$$\begin{aligned} \phi (\mathbf{x},t)= & {} \int _{\Omega } G(\mathbf{x},\mathbf{y}) \, \nabla _y G_D(\mathbf{y},\mathbf{x}_0,t) \cdot \widetilde{\mathbf{f}}_0 \, d \mathbf{y} \nonumber \\+ & {} \int _{\Omega } G(\mathbf{x},\mathbf{y}) \, d \mathbf{y} \int _0^t d \tau \nonumber \\&\,&\int _{\partial \Omega } \psi (\varvec{\xi },\tau ) \, \nabla _y E_t(\mathbf{y}, \varvec{\xi },t-\tau ) \cdot \mathbf{n}_e(\varvec{\xi }) \, d S(\varvec{\xi }) \end{aligned}$$
(98)

and the boundary equation determing \(\psi (\varvec{\xi },t)\) is given by

$$\begin{aligned} \psi (\varvec{\xi },t)= m_t(\varvec{\xi },t) + \int _0^t d \tau \int _{\partial \Omega } K_t(\varvec{\xi },\varvec{\eta },t-\tau ) \, \psi (\varvec{\eta },\tau ) \, d S(\varvec{\eta }) \end{aligned}$$
(99)

for \(\varvec{\xi } \in \partial \Omega\), where

$$\begin{aligned} m_t(\varvec{\xi },t)= \int _{\Omega } \mathbf{n}_e(\varvec{\xi }) \cdot \nabla _x G(\mathbf{x},\mathbf{y}) |_{\mathbf{x}=\varvec{\xi }} \, \nabla _y G_D(\mathbf{y}, \mathbf{x}_0,t) \cdot \widetilde{\mathbf{f}}_0 \, d \mathbf{y} \end{aligned}$$
(100)

and

$$\begin{aligned} K_t(\varvec{\xi },\varvec{\eta },t)= & {} \int _{\Omega } \mathbf{n}_e(\varvec{\xi }) \cdot \nabla _x G(\mathbf{x},\mathbf{y}) |_{\mathbf{x}=\varvec{\xi }} \nonumber \\&\,&\nabla _y E_t(\mathbf{y}, \varvec{\eta },t) \cdot \mathbf{n}_e(\varvec{\eta }) d \mathbf{y} \end{aligned}$$
(101)

6 Concluding remarks

In this article we have presented a physical approach towards the determination of the hydrodynamic Green functions for steady and unsteady Stokes problems.

While in the case of steady flows, this analysis provides just “another” derivation of the Oseen tensor, the application to unsteady Stokes problems yields a simple and “physically readable” expression for the corresponding Green functions and, more importantly, new insights on the physics of the hydrodynamic propagation.

The gradient gauge decomposition, representing the starting point of the present mathematical analysis of the Green functions, is the simplest way to decouple the various agents determining the space-time properties in the solution of the unsteady Stokes equations. The latter can be always viewed as the superposition of a vector field \(\mathbf{v}^\prime (\mathbf{x},t)\), the properties of which are those of any scalar field evolving according to a parabolic equation in space and time, and of a gradient gauge, that instantaneously enforces the condition of vanishing divergence. Even for unsteady Stokes flows, the equation for the gradient gauge \(\phi (\mathbf{x},t)\) is timeless, leading to a Poisson equation, the fundamental properties of which determine at any time \(t>0\) a non-vanishing velocity field \(\mathbf{v}(\mathbf{x},t)\) that decays as a power law of the distance from the point of application of the initial impulsive forcing. This phenomenon can be simply interpreted by resolving a vectorial impulsive forcing into its two, irrotational and solenoidal, Helmholtz components.

The power-law decay of the gradient gauge determines a physical inconsistency in the long-term properties of the velocity field in incompressible unsteady Stokes flow, and the paradoxical phenomenon of infinite propagation velocity of viscous stresses, the origin of which is twofold: (i) the Newtonian relation connecting the stress and the deformation tensors, leading to the parabolic evolution equation for \(\mathbf{v}^\prime (\mathbf{x},t)\), and (ii) the incompressibility condition which, independently of the evolution equation adopted for \(\mathbf{v}^\prime (\mathbf{x},t)\), generates the unphysical large-distance power-law scaling (48) for the velocity field \(\mathbf{v}(\mathbf{x},t)\).

In the resolution of this paradox, the incompressibility condition represents the unphysical constraint determining non-locality in the viscous stress propagation. In point of fact, the resolution of this paradox is not only of theoretical interest but admits many important applications related to hydrodynamics at small lengthscales. The rapid progress in microfluidics and the detailed experimental and theoretical investigations on the motion of colloidal particles at microscale (Brownian motion) [32, 33] push the hydrodynamic theory towards a more physically consistent description of fluid-particle interactions, overcoming the manifestly paradoxical inconsistency in the large-scale behavior of the velocity fields. In this framework, the generalization of the evolution equation for the velocity field of a liquid phase, possessing dissipative viscous propagation of the internal stresses, and a finite propagation velocity for the density and pressure fields, leading as a byproducing to a stress propagation with bounded velocity is the missing link in the formulation of a low-Reynolds number hydrodynamics of “almost incompressible” flows, in which incompressibility emerges as a steady-state property once dynamic wave-like motion in the liquid has been set down. A formulation of the Navier–Stokes equations resolving these problems can be found in [36]. This problem can also be tackled via experimental analysis, by considering e.g. the perturbation induced by the motion of a micrometric particle, trapped in an optical trap or following a prescribed trajectory using optical tweezers and analyzing via \(\mu\)-PIV equipment the resulting velocity profile in the intermediate and far fields.