Excess Diffusion of a Driven Colloidal Particle in a Convection Array

Excess diffusion is now a well-understood effect,^[1–3] occurring whenever a Brownian particle is transported through an array of geometric constrictions directed along a set direction, e.g., the x-axis. The most investigated case-study is represented by the tilted washboard potential [Fig. 1(a)], whereby increasing the intensity of the longitudinal force, F, close to a threshold value, F_th, causes a particle's sudden depinning. The locked-to-running transition is signalled by a large peak of the particle's diffusion constant, D, in the force direction. Upon lowering the strength of the thermal noise, D₀, depinning occurs as F approaches the maximum confining force exerted by the washboard potential,^[4] F_th L/2π = V₀, and the relative diffusion peak height, D/D₀, diverges proportionally to (V₀/D₀)^−2/3.^[3]

Zoom In Zoom Out Reset image size

Another interesting example of excess diffusion has been detected in the so-called entropic channels.^[5–7] A particle driven along a narrow corrugated channel squeezes its way through bottlenecks (pores), which oppose the action of the external drive with effective barriers of an entropic nature.^[6] In a mirror-symmetrical, two-dimensional (2D) channel with sinusoidal boundaries, ± w(x) [Fig. 1(b)], the particle's diffusion along the channel's axis is closely approximated to a one-dimensional (1D) Brownian motion with "entropic" potential, A(x) = −D₀ ln [2w(x)/Δ],^[6] and effective local diffusion constant, D(x) = D₀/[1 + w'(x)²]^1/3.^[8] This approximation holds good for ε ≪ Δ ≪ L, where ε and Δ denote the pore and channel widths, respectively, and L is the length of the w(x) unit cell. In the presence of a longitudinal force, the above Fick–Jacobs–Zwanzig approximation scheme is tenable, as long as FL/2π ≪ D₀. Nevertheless, numerical simulations show that an excess diffusion peak can emerge for FL/2π ∼ D₀, in cases where the drive is strong enough to overcome the entropic barriers.^[7,9] The main difference between washboard potential and the entropic channel models lies in the nature of the barriers opposing particle transport, i.e., energetic and entropic barriers, respectively.^[10,11] Accordingly, the response to an external bias is of the Arrhenius type for the former, and is non-exponential for the latter. Both classes of transport models are widely studied and utilized in the current literature.^[12–14]

An apparently unrelated problem is Brownian diffusion in convection flows. This is a recurrent problem in today's nanotechnology^[15,16] with relevance to applications in chemical engineering and combustion.^[17] Under quite general conditions, the combined action of advection and thermal fluctuations in the suspension fluid accelerates tracers' diffusion, an effect known as advection enhanced diffusion (AED).^[18,19]

Let us consider an overdamped point-like Brownian particle of unit mass, suspended in a 1D array of counter-rotating convection rolls with periodic stream function [Fig. 1(c)]

$$\begin{eqnarray}&&\psi (x,y)=({U}_{0}L/2\pi)\sin (2\pi x/L)\sin (2\pi y/L),\end{eqnarray} \tag{ 1 }$$

where 0 ≤ y ≤ L/2; x and y are the particle's coordinates, L denotes the size of the flow unit cell, U₀ is the maximum advection speed at the rolls' separatrices, Ω_L = 2π U₀/L is the maximum vorticity at their centers, and D_L = U₀ L/2π denotes an intrinsic flow diffusion constant. Perturbed by thermal fluctuations of strength D₀, the particle undergoes normal diffusion, together with with asymptotic diffusion constant, D, along the array. AED takes place at high Péclet numbers, Pe = D_L/D₀ ≫ 1, whereby D is discovered to be larger than the free diffusion constant, D₀. This effect has been explained in Refs. [18–21], where it is noted that, at low noise, an unbiased particle jumps between convection rolls by virtue of advection, which drags it along the rolls' outer flow layers, or flow boundary layers (FBL), centered around the rolls' separatrices. Thermal diffusion across such narrow FBLs favors the particle's roll-jumping, thereby enhancing spatial diffusion.^[22] In view of this mechanism, AED alone cannot be considered to be an excess-diffusion effect, and can be expected to occur in both square^[18,19] and linear convection arrays.^[20,21]

In spite of their common 1D geometry, the three diffusion setups, as shown in Fig. 1, exhibit substantially different transport properties. For instance, the 1D probability density function (pdf), p(x), in the washboard potential is spatially modulated, in accordance with a Boltzmann distribution. In an entropic channel, the 2D pdf, p(x,y), is uniform, whereas the relevant p(x) is not (as would be apparent in a Fick–Jacobs–Zwanzig approximation). Furthermore, their mobility function, μ = $\langle \dot{x}\rangle$/F, as a function of F, increases from 0 to 1, sharply around F_th for the washboard potential, and smoothly for F ∼ 2π D₀/L in a narrow sinusoidal channel, with the excess diffusion peak being more pronounced in the former. Significantly, in a convection array, both p(x,y) and p(x) are uniform, and μ = 1, regardless of the value of F. This raises a question as to whether a longitudinal bias may affect diffusion in a convection array.

Here, we investigate a distinct instance of advection-assisted diffusion, which manifests its own distinct, huge, excess diffusion peak. When applying a longitudinal force of increasing intensity, the vertical flow layers separating the individual convection rolls rapidly become ineffective. In contrast, the flow layers parallel to the array edges begin to play a more prominent role in the diffusive process. Therefore, even if the spatial distribution of the particle in the array is uniform, its diffusion is advection-enhanced along the top and bottom edges of the array. This produces a huge excess diffusion peak for values of force intensity proportional to the maximum advection speed, U₀.

Model. A Brownian particle diffusing in the convection array of streaming function ψ(x,y) of Eq. (1) obeys the Langevin equations,

$$\begin{eqnarray}&&\dot{x}={u}_{x}+F+{\xi }_{x}(t),\dot{y}={u}_{y}+{\xi }_{y}(t),\end{eqnarray} \tag{ 2 }$$

where u = (u_x , u_y ) = (∂_y , − ∂_x )ψ is the incompressible advection velocity vector, ∇ · u = 0, and F is a tunable, uniform longitudinal force. As illustrated in Fig. 1(c), the array unit cell comprises two counter-rotating convection rolls. The random sources, ξ_i (t) and i = x,y, represent stationary, independent, delta-correlated Gaussian noise, ⟨ξ_i (t)ξ_j (0)⟩ = 2 D₀ δ_ij δ (t), modeling equilibrium thermal fluctuations in a homogeneous, isotropic medium. For a particle of unit mass, D₀ coincides with its free diffusion constant in the absence of advection. Since the flow parameters, L and ${\varOmega }_{{\rm{L}}}^{-1}$, can be adopted as convenient length and time units, the only remaining tunable parameters in our model are the noise strength, D₀ (in units of D_L), and the bias, F (in units of U₀). The stochastic differential Eqs. (2) are numerically integrated by means of a standard Mil'shtein scheme^[23] (Fig. 2). To ensure numerical stability, the numerical integrations are performed using a very short time step, 10⁻⁵ − 10⁻⁴. The stochastic averages reported here have been computed over at least 10⁷ samples (trajectories). Computing the asymptotic diffusion constant, D = lim_{t → ∞} ⟨ Δ x²(t)⟩/2t, with Δx(t) = x(t) − ⟨x(t)⟩, requires extra caution, in that, at high Péclet numbers, the advected particle may take an exceedingly long time to exit a convection roll.

Zoom In Zoom Out Reset image size

At zero bias, the asymptotic diffusion constant, D, changes from $D=\kappa \sqrt{{D}_{{\rm{L}}}{D}_{0}}$ for Pe ≫ 1 to D = D₀ for Pe ≪ 1. The constant, κ, depends on the convection array's geometry and boundary conditions.^[18,21] For the unbiased, square free-boundary convection array given in Eq. (1), κ ≃ 1.07,^[18] which is consistent with previous numerical results.^[22] The crossover between these diffusion regimes is clearly localized around D₀ ≃ D_L;^[24] as such, AED occurs where D₀ < D_L. Many numerical and experimental papers support the FBL-based interpretation of AED.^[25–29]

Results. We begin our investigation of Brownian transport in the convection array given by Eq. (1) by qualitatively analyzing the trajectories of a tracer driven by a longitudinal force of increasing intensity, F (Fig. 2). At F = 0, panel (a), the FBL structure is apparent. The tracer is advected along the rolls' boundaries; after completing several rounds, it jumps into an adjacent roll. Occasionally, it gets stuck in the so-called stagnation corners, i.e., around those points where the flow velocity, u , vanishes. The FBL width, δ, is of the order of the length diffused by a free Brownian particle during a full convection round, where δ = (D₀/Ω_L)^1/2.^[18,19,22]

On increasing F, as shown in panel (b), the tracer's circulation along the rolls' boundaries ceases, and its retrapping inside a convection cell becomes infrequent. Moreover, the transition between the dynamical regimes of panels (a) and (b) occurs in the vicinity of a certain value, F_c, of the drive. We estimated F_c by equating the FBL width, δ, to the net longitudinal displacement undergone by the driven tracer as it crosses the array, F/4Ω_L, so that

$$\begin{eqnarray}&&{F}_{{\rm{c}}}/{U}_{0}=4\sqrt{{D}_{0}/{D}_{{\rm{L}}}}.\end{eqnarray} \tag{ 3 }$$

Further numerical evidence of the critical nature of the dynamical transition taking place at F ∼ F_c will be provided in an upcoming full-length paper. Here, we simply remark that the FBL breaking mechanism induces oscillating trajectories of the highest amplitude. Transport across the convection rolls is governed primarily by advection, with bias merely serving to kick the tracer across an FBL, and into the next roll. Retrappings and multiple rounds of the tracer inside a roll are required to ensure that the average drag speed is $\langle \dot{x}\rangle$ = F, i.e., μ = 1, for any value of F.

On further increasing F, as illustrated in panel (c), the drive eventually becomes equal to the advection drag. Under these conditions, retrapping of the tracer becomes very rare. The corresponding drive value, F_m, can be estimated by comparing the bias drag time across a ψ(x,y) unit cell, where τ_F = L/F, to the advection drag time along the largest circle of radius L/4 enclosed in a roll, such that τ_ψ = (π L/2)/U₀; as such,

$$\begin{eqnarray}&&{F}_{{\rm{m}}}/{U}_{0}=(2/\pi).\end{eqnarray} \tag{ 4 }$$

Finally, for even larger drive values, as depicted in panel (d), the driven particle crosses each roll so fast that advection has little time to perturb its trajectory, the effect of advection is averaged out, and the trajectories begin to resemble those of a driven Brownian particle at zero flow.

This qualitative description is supported by the curves of the autocorrelation function, C(t) = ⟨y(t)y(0)⟩, plotted in Fig. 3 for increasing drive intensities. The limiting values for all curves are as follows: C(0) = ⟨y²⟩ = L²/12 and C(∞) = ⟨y⟩² = L²/16. All exhibit oscillating transients with angular frequency, ω(F), which, as shown in the inset, appears to increase with F, in accordance with the linear law,

$$\begin{eqnarray}&&\omega (F)=\omega (0)+{\varOmega }_{{\rm{L}}}(F/{U}_{0}).\end{eqnarray} \tag{ 5 }$$

Here, ω(0) is the average circulation frequency inside a convection roll. Its numerical estimate, ω(0) = 0.359, closely agrees with an earlier analytical estimate in the noiseless limit.^[30] For exceedingly large values of F, F ≫ U₀, ω(F) ≃ 2π F/L, as would be expected for straight trajectories "hitting" the vertical roll separatrices over a time period τ_F . Based on the fitting law (5), the optimal synchronization of bias and advection time modulations occurs for $F={{F}{^{\prime} }}_{{\rm{m}}}$, with $\omega ({{F}{^{\prime} }}_{{\rm{m}}})={\varOmega }_{{\rm{L}}}$, i.e., ${{F}{^{\prime} }}_{{\rm{m}}}/{U}_{0}=1-\omega (0)/{\varOmega }_{{\rm{L}}}\simeq 0.641$, which is comparable with F_m ≃ 0.637 from Eq. (4).

Zoom In Zoom Out Reset image size

The C(t) curves of Fig. 3 suggest that that the oscillatory transients are dominated by the roll-jumping mechanism. The relevant trajectories run along the faster top and bottom FBL branches, as shown in Figs. 2(b) and 2(c), where advection and bias drag are parallel. At a given distance from the array bottom, y, the total longitudinal drag velocity is u_x (y) + F [first Eq. (2)]; averaging the period of the trajectory oscillations over y yields the approximate formula of Eq. (5).^[30]

Another important feature emerging from Fig. 3 is that the amplitude of the transient oscillations reaches its maximum in correspondence with the FBL breakup; i.e., for F ∼ F_c, Eq. (3) yields, for the simulation parameters used here, F_c ≃ 0.126, which is consistent with the data plotted in the figure. In fact, for F ∼ F_c, the particle is being advected along the roll boundaries for the most part, regardless of whether it is trapped in one roll, or jumps between two adjacent rolls.

However, advection also continues to play a role in the diffusion process, even at relatively large F values. At the array's edges, the particle is subjected to a longitudinal advection drag, modeled by the cosine potential ± D_L cos(2π x/F), see first Eq. (2) for y = 0 or π. Given that edge advection acts parallel or antiparallel to the drive, alternately, one would expect a substantial enhancement of the tracer's longitudinal diffusion along the array's edges. On the other hand, the relevant PDFs, p(x,y) and p(x), are both uniform, and the mobility resembles μ = 1 at zero bias. This is because the particle is not pinned against the array's edges, but can diffuse from one edge to the other. Therefore, any enhancement of its diffusion constant should not be mistaken for the excess-diffusion effect observed in a tilted washboard potential [Fig. 1(a)].

The role of advection in the longitudinal diffusion process is illustrated in Fig. 4, where we plot D versus F, for different values of D₀ (at high Péclet numbers). On closer inspection, a few remarkable properties become apparent:

• For F/U₀ → 0 and F/U₀ → ∞, the diffusion curves approach free-diffusion asymptotes: $D(0)=\kappa \sqrt{{D}_{0}{D}_{{\rm{L}}}}$, and D(∞) = D₀, respectively. The former limit corresponds to the well-established AED phenomenon. The latter limit sets in when the bias is so strong as to supersede advection. As anticipated above [Fig. 2(d)], the tracer can then be assimilated to a driven free Brownian particle. Therefore, its diffusion constant tends to D₀.^[4]
• The curves of D versus F peak around the value F^*, which is apparently insensitive to noise level, given that D₀ varies by more than two orders of magnitude. Having related excess diffusion to the tracer's diffusion between the top and bottom edges of the array, we would expect the excess diffusion peak to be located at around F_m, where the difference between the net drag speeds at the edges is the greatest.^[31]
• One remarkable feature is the magnitude of the excess-diffusion peak, D_max, which appears to scale similarly to ${D}_{0}^{-1}$, as illustrated in the inset of Fig. 4. The relevant diffusion ratio, D(F_m)/D(0) ∝ (D_L/D₀)^−3/2, is far more pronounced than if the particle were constrained to move along either of the array's edges. If that were the case, its motion would be governed by an effective 1D washboard potential, with D(F_m)/D(0) ∝ (U₀/D₀)^−2/3.^[3] Instead, our results indicate the presence of a distinct excess diffusion mechanism.
• Finally, the peak profile exhibits a peculiar F dependence: its left (right) side rises (decays) slowly, according to a clear-cut power law, proportional to F² (F⁻²). A working fitting formula for the low-noise excess diffusion peaks of Fig. 4 may be

  $$\begin{eqnarray}&&\displaystyle \frac{D}{{D}_{{\rm{L}}}}=K(\displaystyle \frac{{D}_{{\rm{L}}}}{{D}_{0}}){[\displaystyle \frac{F/{F}^{* }}{1+{(F/{F}^{* })}^{2}}]}^{2},\end{eqnarray} \tag{ 6 }$$

  where K and F^* are fitting parameters. Our best fits in Fig. 4 center the excess diffusion peak at around F^* = 0.60 ± 0.01, close to our estimate of F_m in Eq. (4).

Zoom In Zoom Out Reset image size

The fitting formula in Eq. (6) can be interpreted by means of a heuristic argument. The advection flows along the top and bottom edges of the array are spatially modulated, with a period L, and are opposite in phase. Moreover, the subtracted particle's speed, ${\bar{u}}_{x}=\langle \dot{x}\rangle -F$, is modulated by the combined action of advection and drive, respectively, with time constants of τ_ψ , and τ_F [Fig. 2(c)]. Following the analytical approach employed in Refs.[32,33], one can predict that ${\bar{u}}_{x}/{U}_{0}\propto {\tau }_{\psi }{\tau }_{F}/({\tau }_{\psi }^{2}+{\tau }_{F}^{2})$. The particle is not pinned specifically to either edge of the array; rather, it switches between them by diffusing across the array, with the time constant τ_y = (L/2)²/2D₀. This means that the edge switching contribution to D shall be of the order of ${\tau }_{y}{\bar{u}}_{x}^{2}/2$, which is consistent with the fitting formula of Eq. (6).

Finally, as implied in the above discussion, and confirmed via numerical simulation, the reported excess diffusion also occurs along the main axes of the square convection array of Eq. (1), as long as the drive is applied parallel to them.

The main conclusions of this study can be summarized as follows: we find that the FBL mechanism responsible for AED in convection rolls at high Péclet numbers is quite weak in terms of fluid mechanics. In the presence of biases with a modulus above a critical value (proportional to the square root of the noise strength), the FBLs break up. However, by applying a longitudinal bias of modulus comparable with the advection speed, we were able to obtain a huge excess diffusion peak. The numerical evidence presented here refers to an example involving free-boundary convection arrays. Further results, to be published in a forthcoming paper, will show that this conclusion also holds good for rigid-boundary arrays.