Onset of anomalous diffusion in colloids confined to quasi-monolayers

TL;DR

This paper investigates how colloidal particles transition from normal to anomalous diffusion when confined in quasi-monolayers, considering the effects of out-of-plane fluctuations and hydrodynamic interactions.

Contribution

It provides a theoretical and numerical analysis of the crossover from normal to anomalous diffusion in quasi-monolayers with out-of-plane particle fluctuations.

Findings

Identification of the conditions leading to anomalous diffusion

Quantitative description of the crossover behavior

Impact of out-of-plane fluctuations on collective diffusion

Abstract

It has been recently shown that a colloidal monolayer, e.g., formed at a fluid interface or by means of a suitable confining potential, exhibits anomalous collective diffusion. This is a consequence of the hydrodynamic interactions mediated by the three-dimensional (3D) ambient fluid when the particles are confined to reside on a two-dimensional (2D) manifold. We study theoretically and with numerical simulations the crossover from normal to anomalous diffusion as the particles are, in real systems, confined by a 3D external potential and thus have the possibility to fluctuate out of the 2D manifold, thus forming actually a quasi-monolayer.

Equations

$$\frac{D(k)}{D_0}=\sum _{n=-\infty }^{\infty }{\beta }_n{k}^n,$$

$${\mathbf{v}}_i=\sum _{j=1}^N{\mathcal{M}}_{ij}\cdot {\mathbf{f}}_j.$$

$$\frac{\partial P}{\partial t}=\sum _{ij}{\nabla }_{{\mathbf{x}}_i}\cdot \left({\mathcal{M}}_{ij}\cdot {\mathbf{\Phi }}_j\right),$$

$${\mathbf{\Phi }}_j(\{\mathbf{x}\}):=k_{\mathrm{B}}T{\nabla }_{{\mathbf{x}}_j}P(\{\mathbf{x}\})+P(\{\mathbf{x}\}){\nabla }_{{\mathbf{x}}_j}U(\{\mathbf{x}\}),$$

$$U(\{\mathbf{x}\})=U^{\mathrm{int}}(\{\mathbf{x}\})+\sum _{k=1}^{N}V({\mathbf{x}}_k),$$

$${\dot{\mathbf{x}}}_i=\sum _{j=1}^N\left[-{\mathcal{M}}_{ij}\cdot {\nabla }_{{\mathbf{x}}_j}U+k_{\mathrm{B}}T{\nabla }_{{\mathbf{x}}_j}\cdot {\mathcal{M}}_{ij}\right]+{\mathbf{\eta }}_i,$$

$$⟨{\eta }_{i,\alpha }(\{\mathbf{x}\},t){\eta }_{j,\beta }(\{\mathbf{x}\},t^{\prime })⟩=2k_{\mathrm{B}}T{\left[{\mathcal{M}}_{ij}(\{\mathbf{x}\})\right]}_{\alpha \beta }\delta (t-t^{\prime })$$

$$\rho (\mathbf{x},t)=\int d^3{\mathbf{x}}_2\dots d^3{\mathbf{x}}_{N}P({\mathbf{x}}_1=\mathbf{x},{\mathbf{x}}_2,\dots {\mathbf{x}}_N,t).$$

$${\mathcal{M}}_{ij}(\{\mathbf{x}\})=\Gamma {\delta }_{ij}\mathcal{I}+\Gamma (1-{\delta }_{ij})\mathbf{\omega }({\mathbf{x}}_i-{\mathbf{x}}_j),$$

$$\Gamma =\frac{1}{3\pi \eta {\sigma }_{\mathrm{H}}}$$

$$\mathbf{\omega }(\mathbf{x})=\frac{3}{8}\frac{{\sigma }_{\mathrm{H}}}{x}\left(\mathcal{I}+\frac{\mathbf{xx}}{x^2}\right).$$

$$\frac{\partial \rho }{\partial t}=D_0{\nabla }_{\mathbf{x}}^2\rho -{\nabla }_{\mathbf{x}}\cdot \left[\rho \left(\mathbf{u}-\Gamma {\nabla }_{\mathbf{x}}V\right)\right],$$

$$D_0=\Gamma k_{\mathrm{B}}T=\frac{k_{\mathrm{B}}T}{3\pi {\sigma }_{\mathrm{H}}\eta },$$

$$\mathbf{u}(\mathbf{x})$$

$$=-\Gamma \int d^3{\mathbf{x}}^{\prime }\rho ({\mathbf{x}}^{\prime })[{\nabla }_{{\mathbf{x}}^{\prime }}V({\mathbf{x}}^{\prime })]\cdot \mathbf{\omega }(\mathbf{x}-{\mathbf{x}}^{\prime }),$$

$$V=k_{\mathrm{B}}T{\left(\frac{z}{{\ell }_{\text{c}}}\right)}^2,$$

$$\mathbf{x}=\mathbf{r}+z{\mathbf{e}}_z,$$

$$P(\{\mathbf{x}\})=P^{(2D)}(\{\mathbf{r}\})\prod _{i=1}^N\delta (z_i),$$

$$\frac{\partial P^{(2D)}}{\partial t}=\sum _{ij}{\nabla }_{{\mathbf{r}}_i}\cdot \left\{{\mathcal{M}}_{ij}\cdot \left[k_{\mathrm{B}}T{\nabla }_{{\mathbf{r}}_j}{P}^{(2D)}+P^{(2D)}{\nabla }_{{\mathbf{r}}_j}{U}^{\mathrm{int}}\right]\right\},$$

$${\dot{\mathbf{r}}}_i=\frac{1}{8\pi \eta }\sum _{j=1,j\ne i}^N\frac{{\mathbf{r}}_i-{\mathbf{r}}_j}{|{\mathbf{r}}_i-{\mathbf{r}}_j{|}^3}+{\mathbf{\xi }}_i,$$

$$\begin{gathered} ⟨{\xi }_{i,\alpha }(\{\mathbf{r}\},t){\xi }_{j,\beta }(\{\mathbf{r}\},t^{\prime })⟩=2D_0[{\delta }_{ij}{\delta }_{\alpha \beta } \\ +(1-{\delta }_{ij}){\omega }_{\alpha \beta }({\mathbf{r}}_i-{\mathbf{r}}_j)]\delta (t-t^{\prime }). \end{gathered}$$

$$\frac{\partial {\rho }^{(2D)}}{\partial t}=D_0{\nabla }_{\mathbf{r}}^2{\rho }^{(2D)}-{\nabla }_{\mathbf{r}}\cdot ({\rho }^{(2D)}\mathbf{u}),$$

$$\mathbf{u}(\mathbf{r})$$

$$\frac{D(k)}{D_0}=1+{\rho }_0^{(2D)}\frac{\mathbf{k}}{k}\cdot \mathrm{FT}[\mathbf{\omega }]\cdot \frac{\mathbf{k}}{k}=1+\frac{1}{L_{\mathrm{hydro}}k},$$

$$L_{\mathrm{hydro}}:=\frac{4}{3\pi {\sigma }_{\mathrm{H}}{\rho }_0^{(2D)}}$$

$${\rho }_{\mathrm{eq}}(z)=\frac{{\rho }_0^{(2D)}}{\sqrt{\pi }{\ell }_{\text{c}}}{\mathrm{e}}^{-V(z)/k_{\mathrm{B}}T},$$

$${\rho }_0^{(2D)}={\int }_{-\infty }^{+\infty }dz{\rho }_{\mathrm{eq}}(z).$$

$$\rho (\mathbf{r},z,t)={\rho }_{\mathrm{eq}}(z)\left[1+\epsilon (\mathbf{r},z,t)\right].$$

$$\frac{\partial \epsilon }{\partial t}=D_0\left[{\nabla }^2+\frac{{\partial }^2}{\partial z^2}\right]\epsilon +\frac{\Gamma }{D_0}\frac{dV}{dz}\left[-D_0\frac{\partial \epsilon }{\partial z}+{\mathbf{e}}_z\cdot \mathbf{u}\right],$$

$$\begin{gathered} \mathbf{u}(\mathbf{r},z,t)=-\Gamma \int d^2{\mathbf{r}}^{\prime }{\int }_{-\infty }^{+\infty }d{z}^{\prime }{\rho }_{\mathrm{eq}}(z^{\prime })\frac{dV}{d{z}^{\prime }}(z^{\prime })\epsilon ({\mathbf{r}}^{\prime },z^{\prime },t) \\ \times {\mathbf{e}}_z\cdot \mathbf{\omega }(\mathbf{r}-{\mathbf{r}}^{\prime }+(z-z^{\prime }){\mathbf{e}}_z). \end{gathered}$$

$$\epsilon (\mathbf{r},z,t)=\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}{\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{r}}\sum _{n=0}^{\infty }{H}_n\left(\frac{z}{{\ell }_{\text{c}}}\right)c_n(\mathbf{k},t).$$

$${\int }_{-\infty }^{+\infty }dz\rho (\mathbf{r},z,t)={\rho }_0^{(2D)}\left[1+{\delta }_0(\mathbf{r},t)\right],$$

$${\delta }_0(\mathbf{r},t):=\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}{\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{r}}{c}_0(\mathbf{k},t).$$

$$\frac{{\ell }_{\text{c}}^2}{D_0}\frac{d{c}_n(\mathbf{k},t)}{dt}=-[({\ell }_{\text{c}}k{)}^2+2n]c_n(\mathbf{k},t)+\frac{{\psi }_n(\mathbf{k},t)}{L_{\mathrm{hydro}}k},$$

$${\psi }_n(\mathbf{k},t)=-\frac{1}{\pi n!2^n}\sum _{m=0}^{\infty }{\Omega }_{nm}({\ell }_{\text{c}}k)c_m(\mathbf{k},t).$$

Full text

I Introduction

Particle–laden fluid interfaces are a common subject in soft matter physics, and offer an interesting approach to effectively 2D systems for theory and experiment. In many cases, it is a good simplification to treat the system of fluids, interface and particles as a genuine 2D problem, e.g., in order to explore phase transitions Zahn:1997 ; Wille:2002 or clustering behavior in lower dimensions Bleibel:2011 . However, as soon as hydrodynamic interactions are considered, the full 3D nature of the setup becomes important Bleibel:SMC14 . A colloidal monolayer is an example of the configuration of partial confinement, as termed in Ref. Bleibel:SMC14 , because a part of the components of the system is confined (the particles are restricted to move in a plane), whereas other constituents are not confined (the ambient fluid occupies the adjacent volume). The dimensional mismatch between the 2D colloidal dynamics and the 3D hydrodynamic interactions mediated by the ambient fluid flow induces anomalous diffusion (more precisely, superdiffusion) for the collective, i.e., large–scale dynamics of the monolayer NKPB02 ; Bleibel:SMC14 ; Bleibel:JPCM15 ; Dom15 . These theoretical predictions lead to a reinterpretation of experimental results that had actually measured the anomalous collective diffusion Lin:2014 .

The dynamics of the spatial distribution of particles in a colloidal monolayer can be conveniently characterized by a wavenumber dependent diffusion coefficient $D(k)$, that can be expanded in powers of the wavenumber,

[TABLE]

where the constant $D_{0}$ is conventionally taken to be the single–particle diffusion coefficient in the dilute limit. Normal diffusion is then characterized by the absence of negative powers of $k$, so that a large–scale ($k\to 0$) perturbation in the monolayer density would relax $\propto\exp(-k^{2}D_{0}\beta_{0}t)$, exhibiting a Gaussian tail in the distribution of particles in real space. The 3D hydrodynamic interactions lead to a value $\beta_{-1}\neq 0$ in Eq. (1), signaling anomalous diffusion beyond a characteristic length scale $L_{\rm hydro}=\beta_{0}/\beta_{-1}$. This divergence of $D(k\to 0)$ leads, in the real–space particle distribution, to an algebraic decay $\propto x^{-3}$ with the distance $x$ from a density perturbation.

These conclusions rely on constraining the particles to a monolayer, identified conventionally with the plane $z=0$. This is the simplest model of an actual experimental configuration, where the particles are trapped at a fluid interface by wetting forces or are forced to stay within a plane by the effect of a strong external potential (e.g., a sheet–like trap created by optical tweezers, the gravitational field if the particles are sufficiently heavy to reside in a bottom layer Zahn:1997 , or the electrostatic attraction to an interface LvHS07 ). The goal of this work is to relax the strong–confinement assumption by considering the quasi–monolayer configuration. We will allow for a more realistic model in which the position of the particles can fluctuate in the $z$–direction, so that the associated dynamics is truly 3D, and will address how the large–scale anomalous diffusion emerges from the underlying 3D normal diffusion. If the quasi–monolayer is characterized by a small thickness $\ell_{\textrm{c}}$, then the main result of our analysis is that the diffusion is normal on scales below $\ell_{\textrm{c}}$, regardless of the presence of hydrodynamic interactions, but the anomalous–diffusion scenario is observed on scales above $\ell_{\textrm{c}}$.

The article is arranged as follows: in Sec. II the theoretical model is described, including a brief review of the general framework and the emergence of anomalous diffusion. The special case of a harmonic trap in the $z$–direction is studied with detail in the linear approximation in density perturbations. Sec. III presents results from truncated Stokesian Dynamics simulations of particles in the harmonic trap as well as from numerical solutions of the corresponding dynamical evolution equation beyond the linear approximation. The last Section summarizes our conclusions. The Appendices collect the more technical parts of the work.

II Theoretical model

For the theoretical description of the dynamics of colloids, the simplifying assumption can be done that the evolution occurs in the overdamped regime and that the ambient flow can be described with the time–independent Stokes equation, i.e., small Reynolds number and instantaneous adjustment of the flow to the particle configuration. This is usually a good approximation for the diffusive dynamics because the time scale of change of the conserved field “particle density” diverges as the spatial extension of a density perturbations is taken arbitrarily large (but see Ref. Dom15 for a discussion of how the anomalous diffusion is affected by allowing for the dynamical evolution of the ambient flow). Under these approximations, the hydrodynamic interactions mediated by the ambient fluid can be taken completely into account by means of the mobility matrix $\mathcal{M}_{ij}(\{\mathbf{x}\})$, a $3\times 3$ matrix for each particle pair $(i,j)$ Dhon96 ; KiKa91 . It depends on the position $\{\mathbf{x}\}=(\mathbf{x}_{1},\dots,\mathbf{x}_{N})$ of all the $N$ particles forming the colloid and implicitly on their shape and size through the boundary conditions that they impose on the ambient flow. The physical meaning of this matrix is that the velocities $\mathbf{v}_{i}$ of the particles are determined by the forces $\mathbf{f}_{j}$ acting on them as

[TABLE]

Correspondingly, the time evolution of the probability distribution $P(\{\mathbf{x}\},t)$ of a configuration of particles at temperature $T$ is described by the Smoluchowski equation Dhon96 ; Allen:1987 :

[TABLE]

with the potential energy

[TABLE]

consisting of an internal part $U^{\mathrm{int}}$ describing the interparticle forces, and a contribution $V$ by an external single–particle potential (in particular, the potential confining the particles to the plane $z=0$). Equivalently, the dynamical evolution for the individual particle trajectories $\mathbf{x}_{i}(t)$ can be described by the associated Langevin equation Allen:1987 ,

[TABLE]

(the Greek subindices refer to the components of the vectors and tensors.) Particularly relevant is the driving force proportional to the divergence of the mobility: it vanishes in bulk, i.e., when the particles can be distributed in the volume, because the Stokes flow is incompressible. However, in the partial confinement configuration, the particles, but not the fluid are constrained to a plane, and this term has the form of a nonvanishing 2D divergence of a 3D mobility matrix (see, c.f., Eq. (15a)).

The collective (large scale, long time) dynamics of a colloid is described by the evolution of the one-particle density distribution,

[TABLE]

One cannot derive a closed equation for $\rho(\mathbf{x},t)$ from the Smoluchowski equation (3) without the introduction of further approximations because of the multiparticle dependence of both the mobility matrix $\mathcal{M}_{ij}(\{\mathbf{x}\})$ and the potential energy $U(\{\mathbf{x}\})$. The simplest approximation, which will be adopted in this work, is to consider the dilute limit. For the mobility matrix, this implies truncating its expansion at the two–particle level and retaining the asymptotically dominant contributions for large interparticle separations ($\mathcal{I}$ is the identity matrix),

[TABLE]

in terms of the single–particle mobility

[TABLE]

(for our case of spherical particles of diameter $\sigma_{\rm H}$ inmersed in a fluid of viscosity $\eta$), and the Oseen tensor

[TABLE]

Effectively, one is taking into consideration only the longest ranged contribution of the hydrodynamic interactions.

For the potential energy, the dilute limit approximation means $U^{\mathrm{int}}(\{\mathbf{x}\})=0$ in Eq. (4), so that the particles do not interact directly with each other (“ideal gas” approximation). With this approximation and Eq. (7), one can obtain from Eq. (3) the following (nonlinear) evolution equation for the one–particle density:

[TABLE]

after integrating by parts and accounting for the incompressibility constraint $\nabla_{\mathbf{x}}\cdot\bm{\omega}(\mathbf{x})=0$. Since the Oseen tensor is the Green function of the Stokes equation in an unbounded volume, one can interpret the field $\mathbf{u}(\mathbf{x})$ as the ambient flow induced by the external forces acting on the particles. Then, Eq. (10a) describes the evolution of the particle density due to Brownian motion and the simultaneous drag by the ambient flow and the external force.

Although the dilute limit approximation is useful for the purpose of this work, it can be relaxed. Thus, in order to account for the effect of direct interparticle forces, local equilibrium approximations for the potential energy landscape are customary. For instance, the so-called dynamical density functional theory MaTa99 and its extension to include the effect of the hydrodynamic interactions Rex:2009 ; Donev:2014 ; Bleibel:JPCM15 . The relevant result is that, for the large–scale dynamics, the effect of the short–ranged interparticle forces shows up as a (possibly density–dependent) change in the numerical value $D_{0}$ of the diffusion coefficient. Thus, it is not expected that the inclusion of direct interparticle forces will alter the qualitative picture, particularly that concerning anomalous diffusion in the monolayer configuration. This latter expectation is actually confirmed by numerical simulations of monolayers composed of interacting particles (capillary monopoles Bleibel:SMC14 , hard spheres GNK16 , Lennard–Jones particles PPD17 ; see also the discussion in Ref. Bleibel:JPCM15 ).

Similarly, for monolayers formed at or close to a fluid interface, a more realistic description of the mobility matrix is possible that accounts for the different values of the fluid viscosity and the particle positioning off the interface Jones:1975 ; AdBl78 ; LoLa14 . Nevertheless, it turns out that the dominant far–field behavior is given again by the Oseen tensor, but with the viscosity $\eta$ in Eq. (7) replaced by the arithmetic mean of the fluid viscosities111The Oseen tensor, decaying as $1/x$, is associated to a so-called “Stokeslet” KiKa91 . The corrections thereof can be written as combinations of “stresslets” and “rotlets” (decaying as $1/x^{2}$) and higher–order terms AdBl78 ; LoLa14 .. Therefore, no qualitative change in the large–scale behavior is expected either.

II.1 The partial confinement configuration

One can consider the particular case that the single–particle external potential $V$ depends only on the $z$–coordinate and has the proper form to force the confinement of the particles within a sheet about $z=0$ of width $\ell_{\textrm{c}}$, see Fig. 1. A good example is a harmonic potential,

[TABLE]

which will be developed in detail below. Before, however, we note some remarks valid beyond this specific form of the potential. One can address two limiting behaviors. In the limit $\ell_{\textrm{c}}\to\infty$ (absence of confinement), one effectively has $V\to 0$ and the Smoluchowski equation (3) describes the dynamics when the particles can explore the whole 3D volume. In the approximated model described by Eqs. (10), one recovers normal diffusion with the single–particle diffusion coefficient $D_{0}$. One notices that the approximations leading to these equations are too simple to capture any effect by the hydrodynamic interactions that would induce a renormalization of the value of $D_{0}$. These appear when short–distance effects are taking into account, like in the more realistic hard–sphere model (see, e.g., Ref. Mazur ). Alternatively, these corrections could be incorporated, in the context of a large–scale, long–time model for the dynamics of the one-particle density, in the form of an effective density–dependence of the single–particle mobility Nozi87 ; Feld88 ; Noet89 ; Lhui90 , which must then be interpreted as a rheological parameter.

The situation is different, however, in the other limiting case, $\ell_{\textrm{c}}\to 0$, which describes perfect confinement of the particles to a monolayer in the plane $z=0$. We introduce the in–plane coordinate $\mathbf{r}=(x,y)$ such that

[TABLE]

and let $\mathbf{k}$ denote the wavenumber vector for the 2D Fourier transform with respect to the in–plane coordinate $\mathbf{r}$. The probability distribution now has the structure

[TABLE]

in terms of the in–plane projected probability distribution $P^{\mathrm{(2D)}}(\{\mathbf{r}\})$. Furthermore, there is no vertical particle current, i.e., $\mathbf{e}_{z}\cdot\bm{\Phi}_{i}=0$ in the Smoluchowski equation (3), and the particle distribution in the $z$–direction is always in equilibrium regardless of the dynamical state of the in–plane distribution. Therefore, upon integrating the Smoluchowski equation over the $z$–coordinates of the particles, one arrives at

[TABLE]

that is, the Smoluchowski equation for the 2D dynamics in the plane $z=0$, where any reference to the confining potential has disappeared. However, the mobility matrix $\mathcal{M}_{ij}(\{\mathbf{r}\})$ still describes a 3D flow (although evaluated at the plane $z=0$), and this dimensional mismatch leads to the anomalous diffusion. The associated Langevin equation for the 2D trajectories $\mathbf{r}_{i}(t)$ of the particles in the monolayer has the form

[TABLE]

when Eq. (5a) is projected onto the monolayer plane with the approximations (7) and $U^{\mathrm{int}}=0$. The force term in Eq. (15a) follows from the observation that $\nabla_{\mathbf{r}_{j}}\cdot\mathcal{M}_{ij}\neq 0$ since the in-plane component of the 3D ambient flow is not 2D incompressible in general, and leads to a force term proportional to $\nabla_{\mathbf{r}}\cdot\omega(\mathbf{r})$ which is formally identical to a Coulombic repulsion, the ultimate cause of the superdiffusive behavior.

The corresponding equation for the 2D particle density field $\rho^{\mathrm{(2D)}}(\mathbf{r})$ is

[TABLE]

Again, any overt signature of the confining potential has disappeared, and Eq. (16a) describes the dynamical evolution driven by the in–plane Brownian diffusion and the drag by the in–plane component $\mathbf{u}(\mathbf{r})$ of the ambient flow. However, the field $\mathbf{u}(\mathbf{r})$ is induced by the in–plane particle current, which is now a relevant source because, unlike in the derivation of Eq. (10), $\nabla_{\mathbf{r}}\cdot\bm{\omega}(\mathbf{r})\neq 0$. As shown in Refs. Bleibel:SMC14 ; Bleibel:JPCM15 , the linearization of this equation for small perturbations about a homogeneous in–plane density $\rho^{\mathrm{(2D)}}_{0}$ yields a wavenumber dependent diffusion coefficient (see Eq. (1)),

[TABLE]

where $\mathrm{FT}[\bm{\omega}]$ denotes the 2D Fourier transform of the 3D Oseen tensor and

[TABLE]

is a characteristic length scale.

II.2 Harmonic confinement: linear theory

The goal is to investigate the transition from one limiting case to the other, with emphasis on the quasi–monolayer configuration. For this purpose we address in detail the linearized theory for the equation (10) with a harmonic confining potential, see Eq. (11). The equilibrium state is given by the Boltzmann distribution,

[TABLE]

where $\rho_{0}^{\mathrm{(2D)}}$ is the projected 2D number density,

[TABLE]

Any particle distribution can be expressed as

[TABLE]

When the model equation (10) is linearized with respect to the small perturbation $\varepsilon$ one obtains an integro–differential equation for the evolution of the fluctuations (see App. A):

[TABLE]

By introducing the Fourier transform for the in–plane $\mathbf{r}$–dependence and an expansion in Hermite polynomials $H_{n}$ for the vertical $z$–dependence, one can write

[TABLE]

Particularly relevant is the coefficient $c_{0}(\mathbf{k},t)$, that describes the Fourier transform of the vertically integrated density profile, which is the 2D density distribution in the partial confinement limit:

[TABLE]

The linearized equation (22) then becomes a set of linear equations for the coefficients $c_{n}(\mathbf{k},t)$ (see App. A):

[TABLE]

Here, $L_{\rm hydro}$ is given by Eq. (18) and the term $\psi_{n}$ encodes the effect of the (long–ranged part of the) hydrodynamic interactions. The dimensionless coefficients $\Omega_{nm}$ are computed in App. B. Of particular relevance is that they are symmetric under the exchange of the indices $n,m$ and vanish when they have different parity. As a consequence, the equations (25) actually form two uncoupled sets of equations: the set for $c_{n}$, $n$ odd, describes particle distributions that are asymmetric in $z$ and whose evolution is driven both by diffusion and the net force exerted by the confining potential. Therefore, we limit ourselves for simplicity to symmetric perturbations in the following, i.e., $\varepsilon(\mathbf{r},-z,t)=\varepsilon(\mathbf{r},+z,t)$, so that the net external force vanishes and we have to consider only the dynamics of the coefficients $c_{n}$ with $n$ even.

The relatively simple structure of Eqs. (25) describes the diffusive relaxation of the modes on a time scale controlled by the length $\ell_{\textrm{c}}$, and the coupling mediated by the hydrodynamic interactions with a strength controlled by the length scale $L_{\rm hydro}$. Two regimes are particularly interesting:

(i) For the very small in–plane scales, $\ell_{\textrm{c}}k\gg 1$ and $L_{\rm hydro}k\gg 1$, neither the confining potential nor the hydrodynamic interactions affect the dynamical evolution appreciably: 3D normal diffusion is recovered because the modes evolve on a time scale $\sim 1/D_{0}k^{2}$ while the term $\psi_{n}$ is strongly suppressed (in addition to the $1/(L_{\rm hydro}k)$ prefactor, it is $\Omega_{nm}(\ell_{\textrm{c}}k\to\infty)\sim 1/(\ell_{\textrm{c}}k)$, see App. B).

(ii) For in–plane scales much larger than the thickness of the quasi–monolayer, $\ell_{\textrm{c}}k\ll 1$, one can recover the scenario originally studied in Ref. Bleibel:SMC14 , as well as derive the leading correction for a nonvanishing thickness $\ell_{\textrm{c}}$. On the one hand, the mode $c_{0}$, associated to the conserved 2D density distribution, is a slow variable, with a characteristic time scale vanishing when $\ell_{\textrm{c}}k\to 0$. On the other hand, the modes $c_{n\geq 2}$ relax on the fast time scale $\sim 2nD_{0}/\ell_{\textrm{c}}^{2}$, signaling the onset of the Boltzmann distribution in the vertical direction. Therefore, the effective dynamics of $c_{0}$ on the slow time scale can be computed approximately by an adiabatic elimination of the fast modes from its equation: the modes $c_{n\geq 2}$ decay to their stationary value at fixed $c_{0}$ and get “enslaved” to the dynamical evolution of the latter. This procedure is detailed in App. C; the final result is

[TABLE]

where $D(k)$ is derived from an expansion in the small parameter $\ell_{\textrm{c}}k$. Therefore, the hydrodynamic interactions give rise to anomalous diffusion for the large in—plane scales satisfying $L_{\rm hydro}k\ll 1$, in agreement with Eq. (17) for the case $\ell_{\textrm{c}}=0$. In addition, there is a finite renormalization of the diffusion coefficient $D_{0}$ for finite values of $\ell_{\textrm{c}}$, but this effect will be hardly observable: it is quantitatively relevant only when $L_{\rm hydro}$ is of the order of $\ell_{\textrm{c}}$, in which case it will be $L_{\rm hydro}k\sim\ell_{\textrm{c}}k\ll 1$, and the anomalous–diffusion effect dominates anyway.

In conclusion, the crossover from 3D normal diffusion to 2D anomalous diffusion occurs smoothly as one shifts the attention from the smallest to the largest scales, with the two length scales $\ell_{\textrm{c}}$ (width of the confining potential) and $L_{\rm hydro}$ (onset of anomalous diffusion) controlling this transition. Our detailed analysis above for the case $\ell_{\textrm{c}}\ll L_{\rm hydro}$ reveals the following hierarchy of dynamical regimes in wavenumber:

[TABLE]

Alternatively, one gets the following scenario in terms of the length scale $r$ of observation: at the smallest scales ($r\ll\ell_{\textrm{c}}\ll L_{\rm hydro}$) the particle distribution diffuses normally inside the quasi–monolayer (regime I); when observed at the intermediate scales ($\ell_{\textrm{c}}\ll r\ll L_{\rm hydro}$), the particle distribution already appears as a perfect monolayer and diffuses normally in the monolayer plane (regime II); and at the largest scales ($\ell_{\textrm{c}}\ll L_{\rm hydro}\ll r$), the monolayer diffusion is anomalously fast (regime III).

III Numerical calculations

III.1 Setup and numerical methods

To further illustrate and complement the results of linear theory for the harmonic confining potential, we performed truncated Stokesian Dynamics (tSD) simulations on the one hand, and solved numerically the density evolution equation (DEE) (10) on the other hand. In both cases the particles were modeled as spheres of diameter $\sigma_{H}=20\;\mu\mathrm{m}$ and the fluid was taken at room temperature ($T=25^{o}\mathrm{C}$) with the viscosity of water $\eta=10^{-3}\;\mathrm{N\,s/m}^{2}$.

The initial particle distribution was constructed as the superposition

[TABLE]

where the background density $\rho_{\mathrm{eq}}(z)$ is given by Eq. (19), which is modified by the radially symmetric overdensity $\Delta\rho(r,z)$. For the latter, we investigated two cases:

• (i)

A planar overdensity which is equilibrated in the $z$–direction but is constant and nonzero on a disk of radius $R$ in the $x$–$y$–plane,

[TABLE]

With this setup we will exemplify the behavior in regimes II and III defined above.

• (ii)

A narrow and isotropic peak of width $l_{\rm G}=\ell_{\textrm{c}}/10$ which is centered at $(r,z)=(0,0)$,

[TABLE]

with the choice $\rho_{\rm G}=(A_{0}/\sqrt{\pi})(\rho_{0}^{\mathrm{(2D)}}R^{2}/l_{\rm G}^{3})$, so that the number of particles in the overdensity is the same as in the planar one (28). This case will allow us to address the behavior in the regime I defined above.

For case (i) we have obtained results from both tSD simulations and the solution of the DEE, while case (ii) has been investigated with the DEE only.

The tSD simulations solve the evolution of a collection of $N$ particles whose dynamics is given by the Langevin equations (5). No direct interaction is considered, $U^{\mathrm{int}}=0$, and the mobility matrix is approximated as in Eq. (7) with the pairwise hydrodynamic interaction given by the Rotne–Prager–Yamakawa tensor RPY:6970 ,

[TABLE]

This tensor is regular at $x=0$ and positive definite, and therefore better suited for particle–based simulations than the Oseen tensor $\bm{\omega}(\mathbf{x})$ (see Eq. (9)). Beyond the dilute limit, a system of particles which do not interact directly (i.e., ideal–gas behavior) but do it hydrodynamically (i.e., with a nonvanishing hydrodynamic radius $\sigma_{\rm H}$ in Eq. (33)) can be realized physically by means of “hairy” particles consisting of a small solid core and a broad polymeric shell around it Lin:1995 ; Bleibel:JPCM15 .

The tSD simulations consisted of $N=1036$ spherical particles. They were performed in a simulation box of extension $L=2000$ $\mu\rm{m}$ in the $x$–$y$–plane with periodic boundary conditions, while unbounded in the $z$–direction since the confinement by the harmonic potential effectively imposes a vanishing particle current at infinity. The initial overdensity in Eq. (28) was simulated with $N_{\rm disk}=188$ particles distributed uniformly within a circular patch of radius $R=100$ $\mu$m and according to the equilibrium profile $\rho_{\mathrm{eq}}(z)$ in the vertical direction. The remaining $N_{\rm b}=848$ particles were distributed similarly but over the whole planar extension of the system. Thus, the effective 2D background density is $\rho_{0}^{\mathrm{(2D)}}=N_{\mathrm{b}}/L^{2}=2.12\times 10^{-4}$ $\mu\mathrm{m}^{-2}$, corresponding to a packing fraction well in the dilute limit, $(\pi/4)\sigma_{H}^{2}\rho_{0}^{\mathrm{(2D)}}\approx 0.07$. This gives $A_{0}=N_{\mathrm{disk}}/(\pi R^{2}\rho_{0}^{\mathrm{(2D)}})=28.2$ in Eq. (28), which represents a large perturbation presumably beyond the scope of the linearized theory. The characteristic length scale of anomalous diffusion (18) associated to this initial configuration was $L_{\rm hydro}=100$ $\mu$m.

The density evolution equation (DEE) (10) is a nonlinear integro–differential equation owing to the hydrodynamic term. We have solved it using an Euler forward scheme for the time evolution. The right hand side of Eq. (10a) was evaluated with Fourier transforms. In the $x$–$y$–plane, radial symmetry was assumed and the corresponding Fourier transforms could be evaluated on an equidistant grid for $\ln r$ using Fast Hankel Transforms. In $z$–direction, we used a Fast Fourier Transform on an equidistant grid. The box size in $z$–direction was chosen $L_{z}=40.96\;R$. Since the system is periodic in $z$–direction through the use of the Fast Fourier Transform, the results for density profiles in radial direction in the $x$–$y$–plane are affected by periodic images for $r\gtrsim L_{z}$.

For later use, we have expressed the thickness of the quasi–monolayer $\ell_{\textrm{c}}$ in terms of a dimensionless confinement parameter

[TABLE]

III.2 Results: diffusion of a planar overdensity

The setup for the planar overdensity given by Eqs. (27, 28) is shown in Fig. 2. This case is a straightforward extension of the planar overdensity investigated earlier in strict 2D confinement Bleibel:SMC14 and focuses on the effect of the finite width of the confining potential in $z$–direction upon the dynamics in the $x$–$y$–plane.

We investigated the range of values $1.25\leq\alpha\leq 10$ for the confinement parameter, corresponding to widths $\ell_{\textrm{c}}$ between 20 $\mu$m and 140 $\mu$m, i.e. for the smallest width the system is close to a monolayer and for the largest width the $z$–extension of the initial overdensity is about as large as the extension in the plane. Of basic interest is the time evolution of the $z$–averaged relative overdensity $\delta_{0}(r,t)$, which corresponds to the inverse Fourier transform of the mode $c_{0}(k,t)$ (see Eqs. (24)). The scales in our setup satisfy $\ell_{\textrm{c}}\lesssim R,\,L_{\rm hydro}$, so it can be conjectured that the expansion of Sec. II.2 in Hermite modes is particularly well suited and fast converging for the density evolution on scales $r\gg R,L_{\rm hydro}$.

In Fig. 3, DEE solutions for $\delta_{0}(r,t)$ are shown for parameters $\alpha=2.5$, 5 and 10 at two different times $t/t_{0}=0.01$ and 0.1 together with the starting configuration. (The time unit $t_{0}=1/(D_{0}\rho_{0}^{\mathrm{(2D)}})$ corresponds to the characteristic Brownian diffusion time at which particles in the plane reach their next neighbor in the background configuration.) As a reference, the relative overdensity in the case of strict 2D confinement (for the same two times, respectively) is shown: The profile decays monotonously in space and shows the instantaneous onset of the $1/r^{3}$–tail characterizing anomalous diffusion. The tail grows in magnitude with time. For a finite thickness of the quasi–monolayer, the spatial density profiles show the same asymptotic, anomalous decay which, however, sets in only at radial distances larger than a critical one. This critical distance also grows with the confinement parameter $\alpha$. This finding has a very straightforward interpretation: Only at distances beyond this critical distance the $z$–confined overdensity appears to be effectively 2D and anomalously decaying. This is in full accordance with the behavior in regime III characterized by the singularity $\propto 1/k$ in the diffusion coefficient, see Eq. (17), derived in the linearized theory.

For radial distances smaller than the critical one we enter regime II. The diffusion of the disklike overdensity at small $r/R<2$ becomes slower with increasing width of the confining potential. At intermediate $r/R\approx 3$ a dip in the overdensity is formed before the profile approaches the anomalous tail for large $r$. This dip is a consequence of the finite thickness of the quasi–monolayer because it is absent in simulations with strict 2D confinement, regardless of the initial extension of the planar overdensity. For the largest width investigated ($\alpha=10$), the overdensity becomes actually negative (i.e., there is a relative depletion) in the dip region. We illustrate this with a time series of overdensity profiles for $\alpha=10$ in Fig. 4 which magnifies the dip region; the overdensity is negative for $1.6<r/R<2.8$.

We observe that the diffusive motion of the edge of the disk becomes slower with increasing $\alpha$; however, it always triggers a hydrodynamic outbound flow at large distances, responsible for the anomalous diffusion and the $1/r^{3}$–tail (regime III). This mechanism drags particles away from the disk edge at a faster rate than they can be replenished by normal diffusion from the disk (regime II), thus developing an initial depletion zone right at the outer edge of the disk, which becomes more conspicuous for larger values of $\alpha$. Interestingly, this means that the 2D effective Green function for the diffusive spread of the overdensity is no longer greater than [math] everywhere (whereas in strict 2D confinement it is).

The initial planar overdensity is equilibrated in $z$–direction and thus only the zeroth Hermite mode $c_{0}(k,t)$ is present. The dynamics, however, leads to a distortion of the Gaussian $z$–dependence and to the appearance of higher Hermite modes $c_{2}(k,t),c_{4}(k,t),...$. These initially grow in time, their strength reaches a maximum at a time $t\sim\ell_{\textrm{c}}^{2}/D_{0}$ and then decays in time. Overall, these higher modes are always much smaller than the leading, zeroth mode.

Next we compare DEE solutions to results from tSD simulations. For a small value $\alpha=1.25$ we show in Fig. 5 DEE solutions and tSD results for the planar overdensity at two times and compare them to the limit of strict 2D confinement and to the limit of 3D Brownian diffusion.

This value of $\alpha$ corresponds to a monolayer of thickness $\approx\sigma_{\rm H}=20$ $\mu$m. We observe that the temporal decay of the overdensity profile is qualitatively as in the strict 2D case but the built-up of the anomalous tail is slower for intermediate $r/R$. At large distances ($r/R>5$), the results from tSD and the DEE solutions (both 2D and 3D) coincide. In both cases, the diffusion of the edge of the disklike overdensity proceeds more slowly as compared to the strict 2D case, but still faster than Brownian Dynamics.

We increase the width to $\alpha=5$ (the width of the confining potential is about $3.5\sigma_{\rm H}=70$ $\mu$m). The tSD profiles clearly confirm the dip which we discussed above for the DEE solutions, see Fig. 6. DEE and tSD profiles agree at large distances, whereas for intermediate $r/R$ the evolution of the profile appears to proceed more slowly in simulations, a fact that may be attributed to the rather small box used in simulations.

However, the main features of the evolution are captured by both methods alike. These main features are the deviations from the case of strict 2D confinement in the depletion zone as well as the onset of anomalous diffusion at larger distances (regime III).

III.3 Results: diffusion of a narrow peak

The setup for the peaklike overdensity given by Eqs. (27, 29) is shown in Fig. 7. The width of the confining potential was set to $\ell_{\textrm{c}}\approx 140$ $\mu$m ($\alpha=10$), so that $l_{\rm G}=\ell_{\textrm{c}}/10=14\;\mu\mathrm{m}$. For $l_{\rm G}\to 0$, the time evolution corresponds to the decay of a $\delta$–peak in the nonlinear DEE (10) (corresponding to a Green function for a linear DEE). This case allows us to study the transition from presumably normal diffusion at small lateral distances (regime I) to anomalous diffusion at larger distances and longer times (regime III). In regime I, we are especially interested in the effect of hydrodynamics in the confined system on smaller length scales; therefore, we have compared the case with hydrodynamic interactions to the 3D, purely Brownian case ($\mathbf{u}=0$ in Eq. (10a)).

Fig. 8(a) shows the time evolution of the overdensity peak in $z$–direction, i.e., $\Delta\rho(0,z,t)$, and Fig. 8(b) shows the evolution in $r$–direction, i.e., $\Delta\rho(r,0,t)$. The initially isotropic Gaussian profile roughly stays Gaussian also at later times but becomes anisotropic. Even though the initial peak is not affected by the confinement, we observe that in the presence of hydrodynamic interactions the peak diffuses faster in lateral $r$–direction than without them (Brownian case), but slower in vertical $z$–direction. This happens already at small times, when the system is still far from being equilibrated in $z$–direction. It can be understood through particle number conservation and the 3D incompressibility constraint ($\nabla_{\mathbf{x}}\cdot\mathbf{u}=0$) that a faster diffusion in $r$–direction must be accompanied with a slower diffusion in $z$–direction: according to Eq. (10), the vertically directed confinement force induces a compressing flow in the $z$–direction and, consequently, an expanding in–plane flow in the $r$–direction.

To elucidate the origin of the faster $r$–diffusion of the blob, we compared the solution to the diffusion of the same peak but without background density. Interestingly, the $r$–diffusion is the same for length scales $r/\ell_{\textrm{c}}\lesssim 1$ (i.e., the background density does not influence it). Only for $r/\ell_{\textrm{c}}\gg 1$ there is a qualitative difference: we observe the anomalous tail in the spatial profile during the decay of the peak on top of the finite background, whereas it is absent in the decay of the peak with no background. As a conclusion, the moderate discrepancy between lateral and vertical diffusion of the peak (an increase of anisotropy) is a hydrodynamic effect in regime I, occurring on scales smaller than the width of confinement.

IV Summary and Conclusions

We have investigated the effect of the hydrodynamic interactions on the collective diffusion of a dilute colloidal suspension when the particles (having a finite hydrodynamic radius) are confined in the vertical direction by a potential of width $\ell_{\textrm{c}}$ in order to build a quasi–monolayer. Hydrodynamic interactions have been approximated by the far–field limit at the two–body level, appropriate for dilute suspensions. The diffusion equation becomes a nonlinear, integro–differential equation in this case. We have investigated collective diffusion using (i) linearized theory, (ii) numerical solutions of the diffusion equation and (iii) truncated Stokesian Dynamics simulations. The analysis of the linearized theory allows the identification of three regimes: On scales much larger than the width of confinement $\ell_{\textrm{c}}$ and the characteristic length $L_{\rm hydro}$, the collective diffusion in the monolayer is always anomalous. For scales below $L_{\rm hydro}$, the density evolution follows 2D normal diffusion, and for scales below $\ell_{\textrm{c}}$, the 3D normal diffusion is recovered. We have confirmed by numerical solutions and simulations that indeed at lateral distances $r$ much larger than the width $\ell_{\textrm{c}}$, the spatial decay of density fluctuations shows instantaneously the signature ($\propto r^{-3}$) of anomalous diffusion characteristic for a perfect 2D monolayer. The numerical approach also allowed the investigation of the transition from 3D to 2D diffusion: at small distances $r\sim\ell_{\textrm{c}}$ one already observes how the hydrodynamic interactions induce faster diffusion in lateral direction but slower in vertical direction. Other peculiar effects also induced by the hydrodynamic interactions, such as the generation of regions of noticeable particle depletion, were observed.

In a recent publication, Panzuela et al. PPD17 address precisely the same problem of the 3D $\to$ 2D crossover in the diffusive dynamics of a monolayer. The numerical simulations presented in Ref. PPD17 are used to measured the in–plane intermediate scattering function of the equilibrium density fluctuations, rather than the decay of macroscopic density profiles, as we have done in this work. Nevertheless, the conclusions agree in both works, which thus represent complementary numerical confirmations of the phenomenology associated to anomalous diffusion in monolayers. For completeness, the detailed relationship between the approach in this work and that in Ref. PPD17 is discussed in App. D.

Our results have far–reaching consequences for the collective diffusion behavior of confined systems in an infinite (or half–infinite medium). These systems encompass bulk colloidal suspensions in an external, sheet–like potential, or colloids and surfactants at fluid interfaces. Whenever lateral distances larger than the confinement width are considered, the collective diffusion must be considered anomalously fast. Experiments on colloidal monolayers indicate a hydrodynamic–induced enhancement of collective diffusion Zahn:1997 . A clear experimental signal for the $1/k$ divergence of the collective diffusion coefficient can be found in Ref. Lin:2014 which is the only experimental work we are aware of. The analysis presented in this work is intended to provide the theoretical framework for the analysis of the 3D $\to$ 2D crossover in possible future experimental realizations of the quasi–monolayer configuration.

It is to be noted that, in many experimental realizations, the relevant configuration is a monolayer in a curved interface, e.g., that of a fluid droplet. Therefore, the extension of the analysis to this case is desirable, with the goal of addressing the effect of curvature on the anomalous diffusion phenomenology. It may be expected that the scenario discussed here is recovered on length scales much smaller than the typical radius of curvature; the general problem is more involved and requires a detailed study.

Acknowledgements.

J.B. thanks the German Research Foundation (DFG) for the financial support through the Project BL 1286/2-1. A.D. acknowledges support by the Spanish Government through Grant FIS2014-53808-P (partially financed by FEDER funds).

Appendix A Linearized equations

Equation (10a) can be rewritten as

[TABLE]

The equilibrium solution given by Eq. (19) implies $\mathbf{u}_{\mathrm{eq}}=0$ (after integrating by parts in Eq. (10)) and

[TABLE]

Therefore, the substitution of Eq. (21) gives

[TABLE]

which leads to the following dynamical equation for $\varepsilon$ after using Eq. (36) and the incompressibility constraint $\nabla_{\mathbf{x}}\cdot\mathbf{u}=0$:

[TABLE]

The linearization of this equation around the unperturbed solution $\varepsilon=0$ follows easily when accounting for the fact that $\mathbf{u}$ is already of linear order in $\varepsilon$ because $\mathbf{u}_{\mathrm{eq}}=0$. In this manner, Eq. (22a) is obtained, while Eq. (22b) is simply the already linear Eq. (10). In particular, for the harmonic confining potential (11) the linearized equation for $\varepsilon$ takes the form

[TABLE]

The expansion (23) can be inverted as

[TABLE]

and Eq. (25a) is obtained from Eq. (38) by using that the Hermite polynomials satisfy

[TABLE]

The function $\psi_{n}(\mathbf{k},t)$ appearing in Eq. (25a) is given as

[TABLE]

which is computed in App. B.

Appendix B Calculation of the coefficients $\Omega_{nm}$

In order to compute $\psi_{n}$ defined by Eq. (40), one first calculates $\mathbf{e}_{z}\cdot\mathbf{u}$ given by Eq. (22b) with the harmonic potential (11) (for the purpose of this Appendix, the explicit time dependence will be dropped from the notation):

[TABLE]

This expression can be evaluated by inserting the 3D Fourier transform of the Oseen tensor KiKa91 ,

[TABLE]

and using that

[TABLE]

so that the integrals over $k_{z}$ and $\mathbf{r}^{\prime}$ can be carried out, which results in

[TABLE]

after using the definition (18) and the expansion (23). Therefore, when this expression is inserted in Eq. (40), one obtains Eq. (25b) with the coefficients

[TABLE]

in terms of the dimensionless quantities

[TABLE]

It is manifest that $\Omega_{nm}(q)$ is symmetric in the indices. Furthermore, a change of variables $\zeta\to-\zeta$, $\zeta^{\prime}\to-\zeta^{\prime}$ in the integrals shows that $\Omega_{nm}(q)$ vanishes if $n$ and $m$ have different parity. It is possible to simplify Eq. (43) and eventually express it in terms of the error function. We are mainly interested, however, in the asymptotic behaviors in the limits $q\gg 1$ and $q\ll 1$, and this can be derived directly from Eq. (43).

When $q\to\infty$, one can evaluate Eq. (43) using Laplace’s formula BeOr78 , because the integral is dominated by the value of the integrand near $\zeta-\zeta^{\prime}=0$. Introducing the new variables $\mu=\zeta+\zeta^{\prime}$, $\sigma=\zeta-\zeta^{\prime}$, one can write

[TABLE]

As $q\to\infty$, this expression can be approximated as

[TABLE]

from where one concludes that $\Omega_{nm}(q)\sim 1/q$.

In the opposite limit $q\to 0$, one can Taylor–expand the integrand in Eq. (43) because $|\mathrm{e}^{-qs}\left(1+qs\right)|\leq 1$ for $s\geq 0$, and thus the integral converges uniformly in $q$. One has

[TABLE]

When this expression is substituted in Eq. (43), the two integrals factorize. They can be computed explicitly by expressing the powers of $\zeta$ and $\zeta^{\prime}$ in terms of the Hermite polynomials and using the associated orthonormality relations:

[TABLE]

For the particular case that both indices $n$, $m$ are even, one obtains

[TABLE]

Appendix C Adiabatic elimination of the fast modes

We introduce the short–hand notations $q:=\ell_{\textrm{c}}k$, $\lambda:=L_{\rm hydro}k$ and $\tau:=D_{0}t/\ell_{\textrm{c}}^{2}$, and define the infinite–dimensional column vector $\mathbf{c}:=(c_{2}\;c_{4}\;\dots)^{\dagger}$, so that the dynamical equations (25) for $n\geq 2$ can be rewritten in compact form as222To avoid a cumbersome notation, we ignore the fact that $n$ and $m$ represent even numbers, but the indices of the components of vectors and matrices must be natural numbers. This should not create ambiguity because the simplicity of the expressions is self-explanatory.

[TABLE]

in terms of the symmetric matrix

[TABLE]

and the column vector $\mathbf{s}:=(\Omega_{20}\;\Omega_{40}\;\dots)^{\dagger}$. The “adiabatic enslaving” of these fast modes (notice that $\mathcal{B}$ does not vanish as $q\to 0$ provided $\lambda\neq 0$) gives the relationship

[TABLE]

Inserting it into the dynamical equation for the slow mode $c_{0}$ (Eqs. (25)) yields

[TABLE]

which becomes the linear Eq. (26a) with the diffusion coefficient

[TABLE]

For consistency with the assumption of “adiabatic enslaving”, this expression is meaningful only in the limit $q\to 0$. From Eqs. (46) and

[TABLE]

[TABLE]

one gets for Eq. (47) the expression

[TABLE]

which is Eq. (26b). The criterion for not retaining higher order terms in this expansion is that they lead to positive powers of $k$ in Eq. (26b). This ultimate goal is also the reason for the careful bookkeeping in powers of both $q$ and $\lambda$ when deriving the expansion.

Appendix D Comparison with Ref. PPD17

Panzuela et al. PPD17 obtain a theoretical result for the short–time collective diffusion coefficient $D^{\mathrm{(short)}}(k)$ by studying the decay of in–plane equilibrium density fluctuations, i.e., the intermediate scattering function

[TABLE]

where $\langle\dots\rangle$ denotes the average over the equilibrium distribution in the initial state, and $\hat{\rho}$ is the microscopic density field. The theoretical analysis assumes a dilute system and focuses onto the short–time regime, i.e., times $t\to 0$ so that the colloidal particles are displaced by an amount much smaller than the mean interparticle separation. In such case, one assumes

[TABLE]

with the definition

[TABLE]

The theoretical model we have devised concerns the time evolution of the average density $\rho=\langle\hat{\rho}\rangle$, see Eq. (6). Nevertheless, the same result can be obtained for the short–time dynamics starting with our linearized equations for the time evolution of a density fluctuation (Eqs. (22)). In terms of the projected average 2D density (see Eqs. (24)),

[TABLE]

one can define the short–time diffusion coefficient as

[TABLE]

Use of Eqs. (25) gives the expression

[TABLE]

This is not a well–defined system–characteristic quantity due to the dependence on the specific initial conditions $c_{m}(\mathbf{k},0)$. However, one can restrict consideration to initial perturbations with $c_{m}=c_{0}\delta_{m,0}$, as is actually done in Eq. (48) when performing the average over equilibrium configurations, for which the different modes $c_{m}$ are uncorrelated: in such case, after evaluating $\Omega_{00}$, see App. B, one obtains

[TABLE]

in terms of the complementary error function $\mathrm{erfc}(q)$. This result coincides exactly with Ref. (PPD17, , Eq. (18)), with the notation $\ell_{\textrm{c}}=\sqrt{2}\,\delta$ (PPD17, , Eq. (1)) and $L_{\rm hydro}=2a/3\phi$ (PPD17, , Eq. (20)). It is to be compared with Eq. (47) derived in the opposite, long–time limit. Both coefficients agree on the dominant, anomalous–diffusion behavior at large scales. Expression (49) is not restricted to the small $q$ limit, but at the price of choosing a certain set of initial conditions ad hoc. Equation (47), on the contrary, is valid only in the limit $q\to 0$, but it incorporates naturally the irrelevance of the initial conditions through the “adiabatic enslaving”.