Active particles under confinement: aggregation at the wall and gradient formation inside a channel

Aggregation of motile organisms, or self-propelled particles in general, is a familiar phenomenon, which is currently driving an intense interest in many scientific disciplines that include physics, biology and sociology¹. Many different physical mechanisms can lead to the aggregation of active systems. For example, these mechanisms can be directional alignment interactions [1–5], density-dependent motility [6, 7] and aggregation at the wall when the system is under confinement [8–11]. In this theoretical paper, I will focus on the latter mechanism: confinement-induced aggregation of active particles. Experimentally, this mechanism is of particular interest given its inevitability due to the finite size of an experimental apparatus. The advent of microfluidic devices as an experimental tool to study active systems further motivates the current study since the channel walls constitute an integral part of the system [12].

In the classical theory of liquids, aggregation of particles close to the wall is well understood and results from the pairwise interactions, such as volume exclusion, of the particles in the system (see e.g. [13] for a review). Density functional theory is an effective way of investigating such a phenomenon, and a parallel effort has been pursued recently for a system of self-propelled rods [9]. On the other hand, if pairwise interactions between the particles become negligible, such a confinement-induced aggregation phenomenon will cease to exist in an equilibrium system, while it will continue to persist in an active system. This is the effect I will focus on here. Specifically, I will investigate a system of active particles that travel at a fixed speed in a channel, undergo rotational diffusion, and react only with the walls. I will derive from first principles a set of equations that govern the density profile of this system, and show that the density function can be written as a series expansion of the Mathieu functions [14]. To make further analytical process, I will introduce a simplified model in which the particles' directions of travel are discretized.

Besides determining quantitatively how confinement induces aggregation, I will also discuss the gradient formation of particle density in a channel that is partitioned into two regions within which the active particles exhibit distinct rotational fluctuations (figure 1(a)). Such a scenario can for instance be engineered by having two background mediums coexisting in a microfluidic channel, and has been recently investigated in a combined experimental and numerical work [15].

Zoom In Zoom Out Reset image size

The structure of this paper is as follows. I will first specify the model under study in section 2, and present three 'predictions' based on intuitive arguments in section 3. The mathematical formulation of the system will be discussed in section 4. A simplified model and its full analysis will be presented in section 5. The paper ends with a discussion of the findings.

I consider a minimal model of active particles in two dimensions, where every particle is assumed to have constant speed u. Noise, of strength a, is incorporated in the direction of travel and the particles interact only with the boundary of the system. Mathematically, the equations of motion for such a N-particle system are

$$\begin{equation} \frac{{\rm d} \mathbf{r}_i}{{\rm d} t} = u \mathbf{v}(\theta_i), \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{{\rm d} \theta_i}{{\rm d} t} = \sqrt{2a}\eta_i(t), \end{equation} \noindent \tag{ 2 }$$

where 1 ⩽ i ⩽ N, R ≡ (r₁,...,r_N), Θ ≡ (θ₁,...,θ_N), v(θ) ≡ (cos θ,sin θ) and (equation (2)) is an Itô stochastic differential equation [16] such that

$$\begin{equation} \langle \eta_i(t) \rangle =0,\quad\langle \eta_i(t)\eta_j(t') \rangle = \delta_{ij} \delta(t-t') . \end{equation} \noindent \tag{ 3 }$$

The specific geometry considered is depicted in figure 1(a) where the walls are assumed to be perfectly rigid and smooth such that the particles are free to rotate at and slide along the walls. Note that in this model, I ignored the thermal translational diffusion D of the particles [17, 18]. As explained in [19], this simplification is valid if D ≪ u²a⁻¹, which applies, e.g. for the bacterium Escherichia coli.

Without loss of generality, I will from now on set the length scale and temporal scale so that the speed u is one and the width of the channel 2d is 2, although I will continue to keep the symbols in the equations so that the dimensions of the terms in the equations are always clear.

I will first analyse the system based on physical arguments in this section. Since there are no interactions between particles by assumption, I can focus on the average behaviour of a single particle. Specifically, denoting the average time the average particle spent at the wall by τ_w, and the average time spend in the channel by τ_c, the boundary condition suggests that if a particle is stuck at a wall, i.e. having its direction of travel pointing towards the wall, then it will only leave the wall if it manages to rotate its direction to point away from it, so the average time spent at the wall can be estimated as the time it takes to rotate, say, by one radian. As a result,

$$\begin{equation} \tau_{\mathrm{w}} \sim 1/a . \end{equation} \tag{ 4 }$$

To estimate τ_c, one has to worry about whether the particles in the channel is dominantly ballistic or diffusive. These two behaviours are dictated by the rotational diffusion constant a. In other words, if a is large, then the particles will have rotated around many times over before reaching the walls and so the motion of the particle is diffusive with an effective diffusion constant D_eff = u²/a (see e.g. chapter 6 in [20]). On the other hand, if a is small, the motion of the particle inside the channel will be predominantly ballistic. Therefore,

$$\begin{equation} \tau_{\mathrm{c}} \sim \left\{ \begin{array}{@{}ll@{}} d^2/D_{\mathrm{eff}}, & {\rm diffusive,}\\ d/u, & {\rm ballistic.} \end{array}\right. \end{equation} \tag{ 5 }$$

Remembering that the units are set such that d = u = 1, τ_c becomes

$$\begin{equation} \tau_{\mathrm{c}} \sim\left\{ \begin{array}{@{}ll@{}} a, & {a \gg 1,}\\ 1, & {a \ll 1.} \end{array}\right. \end{equation} \tag{ 6 }$$

Together with (equation (4)), one may estimate the ratio of the number of particles inside the channel n_c over the number of particles at the walls n_w to be

$$\begin{equation} \frac{n_{\mathrm{c}}}{n_{\mathrm{w}}} \sim \frac{\tau_{\mathrm{c}}}{\tau_{\mathrm{w}}} \sim \left\{ \begin{array}{@{}ll@{}} a^2, & {a \gg 1},\\ a, & {a \ll 1} . \end{array}\right. \end{equation} \tag{ 7 }$$

If the channel is now partitioned into two regions system with distinct rotational diffusion constants a₁ and a₂ (figure 1(a)), an interesting question that is of experimental interests [15] is which side of the region will be of higher density. Focusing on the a₁,a₂ ≫ 1 regime, one may expect that the particles would behave diffusively inside the channel, with effective diffusion constant D₁ ∼ 1/a₁ and D₂ ∼ 1/a₂ in the two distinct regions. Intuitively, if there is a gradient in the system, then the region with the lower diffusion constant would have a higher number of particles since the slow mobility would serve to localize particles [6, 7]. Interestingly, such an effect has also been observed in the embryo of the round worm Caenorhabditis elegans, where regions with lower diffusion constants serve to localize certain diffusing proteins [21, 22]. If this expectation extends to our active system, then one prediction would be that for a₁,a₂ ≫ 1, a₁ < a₂ implies n⁽¹⁾_c < n⁽²⁾_c where n⁽ⁱ⁾_c denotes the number of particles in region i.

Summarizing this section, the following three predictions seem to follow from the intuitive arguments presented above:

• (i)  
  In the symmetric case (a₁ = a₂ = a), if a ≪ 1, then n_c/n_w ∼ a (?)
• (ii)  
  In the symmetric case, if a ≫ 1, then n_c/n_w ∼ a² (?)
• (iii)  
  In the asymmetric case, if a₁ < a₂ and a₁,a₂ ≫ 1, then n⁽¹⁾_c < n⁽²⁾_c (?)

With the mathematical tools developed in the next section, I will show that in fact only the first 'prediction' is true.

I will now provide a mathematically rigorous formulation of the problem. If we denote the probability distribution of the density of particles in the state (R,Θ) at time t by f(t,R,Θ), then the Fokker–Planck equation corresponding to the system is [23]

$$\begin{equation} \frac{\partial f}{\partial t} = \sum_i \left\{a \frac{\partial^2 }{\partial \theta_i^2}f -u\nabla_{\mathbf{r}_i} \cdot [\mathbf{v}(\theta_i) f] \right\}. \end{equation}\noindent \tag{ 8 }$$

Since by assumption, the particles do not interact with each other, the dynamics of the system is fully captured by the corresponding single-particle density function $p(\mathbf {r}_1, \theta _2) \equiv N \int {\mathrm { d}} \mathbf {r}_2 \cdots \mathbf {r}_N\,{\mathrm { d}} \theta _2 \cdots \theta _N f(\mathbf {R}, \Theta )$, which is governed by the following dynamical equation [24, 25]:

$$\begin{equation} \frac{\partial p}{\partial t} = a \frac{\partial^2 p}{\partial \theta^2} -u\nabla_{\mathbf{r}_i} \cdot [\mathbf{v}(\theta)p]. \end{equation}\noindent \tag{ 9 }$$

Focusing on the steady-state solution and assuming the channel is infinite along the y-direction, the equation simplifies to

$$\begin{equation} a\frac{\partial^2 p}{\partial \theta^2} = u\cos \theta \frac{\partial p}{\partial x}. \end{equation}\noindent \tag{ 10 }$$

Employing the separation of variables method by assuming that p(x,θ) = A(θ)B(x), I arrive at the two equations

$$\begin{equation} \partial_{\theta \theta}A = q \cos \theta A , \end{equation}\noindent \tag{ 11 }$$

$$\begin{equation} \partial_x B = qa B, \end{equation} \tag{ 12 }$$

where q corresponds to a set of constants to be determined. Equation (12) can be easily solved

$$\begin{equation}B= K{\rm e}^{Dq x},\end{equation} \tag{ 13 }$$

while equation (11) indicates that the function A is a Mathieu function, usually denoted as ce_2n(θ,q_2n) where n∈N, such that the corresponding characteristic number is zero [14]. Since the characteristic numbers for these types of Mathieu functions are identical for ±q, the expansion for p is

$$\begin{equation} p(x,\theta)= \sum_{n\geqslant 0} \left[a_n\,{\rm e}^{q_{2n}ax} {\rm ce}_{2n} (\theta/2, 2q_{2n}) +b_n{\rm e}^{-q_{2n}ax} \,{\rm ce}_{2n} (\theta/2, -2q_{2n})\right]. \end{equation} \tag{ 14 }$$

As usual, the expansion parameters are obtained by the boundary conditions, which I will determine now.

Denoting the number of particles at the wall on the left (x = −d) travelling in the direction θ by L(θ), L(θ) must be zero for −π/2 < θ < π/2 as particles travelling in those directions would not stay at the wall (figure 2(b)). For π/2 < θ < 3π/2, the equation for L(θ) can be derived by the conservation of particles as follows: the flux of particles that arrive at the wall is −u cos θp^L(θ) where p^L(θ) denotes p(x = −d,θ). At the wall, the particles that are stuck (i.e. particles with direction −π/2 < θ < π/2) undergo rotational diffusion. Therefore, the density function L(θ) satisfies the following differential equation:

$$\begin{equation} a \frac{\partial^2 L}{\partial \theta^2} =-u\cos \theta p^L(\theta) , \quad \pi/2<\theta < 3\pi/2 \end{equation}\noindent \tag{ 15 }$$

with the boundary conditions L(θ = ± π/2) = 0. Finally, the flux of particles that leave the wall corresponds to the flux of particles that reach the angle θ = ± π/2. For instance, for the out-flux source at θ = π/2, the magnitude of the out-flux is |a∂L(θ = π/2)/∂θ|. Together with the out-flux source at θ = −π/2, the overall out-flux at the left wall is

$$\begin{equation} a\left\{\delta(\theta-\pi/2) \left| \frac{\partial L(\theta= \pi/2)}{\partial \theta} \right|+\delta(\theta+\pi/2) \left| \frac{\partial L(\theta=-\pi/2)}{\partial \theta} \right|\right\}. \end{equation} \noindent \tag{ 16 }$$

Zoom In Zoom Out Reset image size

Similarly, the equation governing the steady-state profile at the wall on the right R(θ) is

$$\begin{equation} a \frac{\partial^2 R}{\partial \theta^2} =u\cos \theta p^R(\theta) , \quad -\pi/2<\theta <\pi/2 \end{equation} \tag{ 17 }$$

with the same boundary conditions R(θ = ± π/2) = 0.

As aforementioned, at the steady-state, the particles that are accumulated at the walls eventually re-enter the channel at four delta sources (two at either wall) (see figure 2(a)). To account for these influxes of particles from the walls, equation (10) is modified as follows:

$$\begin{equation} \fl a\frac{\partial^2 p}{\partial \theta^2} = u\cos \theta \frac{\partial p}{\partial x} +a\delta(x-d)\delta(\theta \pm \pi/2) \left| \frac{\partial R}{\partial \theta} \right| +a\delta(x+d)\delta(\theta \pm \pi/2) \left| \frac{\partial L}{\partial \theta} \right|. \end{equation} \noindent \tag{ 18 }$$

In summary, the steady-state density profile of the active particles is obtained by solving the coupled differential equations (15), (17) and (18).

The standard way to proceed at this point is to substitute into the above coupled differential equations the expansion for p in terms of the Mathieu functions (equation (14)), and then try to ascertain the expansion coefficients accordingly. In practice, further analytical process is difficult since manipulations of the Mathieu functions are analytically intractable. Indeed, the numerical challenge to determine the characteristic numbers for Mathieu functions is comparable to performing Brownian dynamics simulation of the system [14] (see footnote ²). Therefore, in order to further reveal the underlying physics of this problem, I will now introduce and analyse a simplified version of the original model that admits full analytical solution.

In this simplified model, the angular components are partitioned into six directions as shown in figure 1(b) and each particle will switch to a neighbouring direction of travel with rate a. I will denote the density function inside the channel travelling towards the ith direction by p_i, and those at the wall on the right and on the left by R_i and L_i respectively. The translational symmetry along the y-direction indicates that p₁ = p₃, p₆ = p₄, and also that the functions p only depend on x but not y. As a result, the dynamical equations for the system are

$$\begin{equation} \partial_t p_1 = a(p_2+p_6-2p_1) -b\partial_x p_1, \end{equation} \tag{ 19 }$$

$$\begin{equation} \partial_t p_2 = 2a(p_1-p_2) -\partial_x p_2, \end{equation} \tag{ 20 }$$

$$\begin{equation} \partial_t p_5 = 2a(p_6-p_5) +\partial_x p_5, \end{equation} \tag{ 21 }$$

$$\begin{equation} \partial_t p_6 = a(p_1+p_5-2p_6) +b\partial_x p_6, \end{equation} \tag{ 22 }$$

where $b=\sqrt {3}/2$ since b = cos π/3. In matrix form, dp/dx = Mp where p = (p₁,p₂,p₅,p₆)^T and

$$\begin{equation} M=\left(\begin{array}{@{}cccc@{}} -2a/b & a/b & 0 & a/b\\ 2a & -2a & 0 & 0\\ 0 & 0 & 2a & -2a\\ -a/b & 0 & -a/b & 2a/b \end{array}\right). \end{equation} \noindent \tag{ 23 }$$

The four eigenvalues associated to the matrix M are {0,0,−aγ/b,aγ/b} where $\gamma = \sqrt {3+4b+4b^2}$, and the corresponding eigenvectors are the column vectors in the following matrix:

$$\begin{equation} V=\left(\begin{array}{@{}cccc@{}} 0 & 1 & 2+\gamma-\frac{2b}{2b+\gamma} & 1\\ 0 & 1 & -4(b+b^2)-\frac{2b(2-\gamma)}{2b+\gamma} & \frac{2b}{2b+\gamma}\\ 0 & 1 & \frac{2b}{2b+\gamma} & -4(b+b^2)-\frac{2b(2-\gamma)}{2b+\gamma}\\ 0 & 1 & 1 & 2+\gamma-\frac{2b}{2b+\gamma} \end{array}\right). \end{equation} \tag{ 24 }$$

Furthermore, the steady-state solution is of the form

$$\begin{equation} p_h(x) = \sum_k c_k\, {\rm e}^{\sigma_k x} V_{hk}, \end{equation} \tag{ 25 }$$

where σ_k∈{0,0,−aγ/b,aγ/b} are the corresponding eigenvalues of M. The above expression is to be compared to equation (14) in the continuum model.

To determine the coefficients c_k in equation (25), I employ the discrete version of the boundary conditions as discussed in the previous section (equation (15)), which for the wall on the left gives

$$\begin{equation} L_1 =L_2 =0, \end{equation} \tag{ 26 }$$

$$\begin{equation} L_5 = \frac{u}{a} \left( p_5^L+bp_6^L\right), \end{equation} \tag{ 27 }$$

$$\begin{equation} L_6 = \frac{u}{a} \left( \frac{p_5^L}{2}+bp_6^L\right) \end{equation} \noindent \tag{ 28 }$$

and the flux of particles out of the wall on the left is (see equation (16))

$$\begin{equation} p^L_1 = \frac{aL_6}{bu} \end{equation} \tag{ 29 }$$

$$\begin{equation} \hspace{1.4pc}= \frac{p_5^L}{2b}+p_6^L, \end{equation} \tag{ 30 }$$

where the second equality follows from equation (28). Similar expressions can be readily found for the R_i. In particular, the condition in equation (30) and the corresponding one for p^R_i enable me to determine the coefficients in the expansion in equation (25) completely.

5.1. Symmetric case (a₁ = a₂ = a)

In the symmetric case, the interesting quantity to investigate is the ratio of the number of particles inside the channel n_c versus the number of particles at the wall n_w. The analytical expression for n_c/n_w can be obtained straightforwardly with the help of Mathematica [26]. Although the analytical expression is complicated and difficult to decipher, it does provide the following asymptotic expressions with respect to a:

$$\begin{equation} \displaystyle \frac{n_{\mathrm{c}}}{n_{\mathrm{w}}} =\left\{\begin{array}{@{}ll@{}} \displaystyle \frac{4}{3 b} \ a \simeq 1.54a , \quad & a \ll 1,\\ \displaystyle \frac{12 a \left(-1+\gamma+2 b \gamma+\gamma^2\right)}{\gamma \left(7+12 b^2+4 \gamma+2 b (7+3 \gamma)\right)} a \simeq 1.17a , \quad & a \gg 1. \end{array}\right. \end{equation} \noindent \tag{ 31 }$$

This result indicates that the ratio increases linearly in both asymptotic regimes, albeit with two different slopes. These analytical results are in perfect agreement with results from Brownian dynamics simulation (figure 2(d)). Referring to section 3, expectation (i) is thus verified while expectation (ii) is invalidated. I will come back to discuss why this is the case in the discussion section.

5.2. Asymmetric case (a₁ ≠ a₂)

In the two-region case as depicted in figure 1(a), the expansion coefficients in equation (25) are obtained by employing the boundary conditions at both walls (equation (30)), and by matching the densities p_i at x = 0. Note that the matching condition at x = 0 ensures the continuity of the density functions but not their slopes. This is fully supported by simulation result (figure 3(c)).

Zoom In Zoom Out Reset image size

Figure 3(c) also shows that there is a gradient of particle density inside the channel, one natural question is which region possesses more particles. The answer can again be obtained from the analytical solution for p_i, R_i and L_i, and the result is that the region with a larger rotational diffusion constant always has a smaller number of particles, whether the comparison concerns the particles inside the channel only (figure 3(a)) or includes the particles at the walls (figure 3(b)). In other words, if a₁ < a₂, then n⁽¹⁾_c > n⁽²⁾_c. Hence, the expectation (iii) discussed in section 3 is not valid. Note that this surprising phenomenon has also been observed previously in a simulation study of a similar type of active system [15].

As we have seen, among the three expectations arisen from intuitive arguments in section 3, only the first expectation is verified. The reason that expectation (ii) is invalid is due to the overall conservation in the number of particles in the system. In other words, although it remains true that as a → ∞, n_w ∼ 1/a, n_c cannot scale with a as this would lead to a divergence in the total number of particles. So as a increase, n_c also increases but eventually saturates to a number bounded below by the total number of particles in the system. This is why n_c/n_w remains linear in a even in the large a regime.

In the asymmetric case, in contradiction to the expectation (iii) from section 3, a higher a₂ in relation to a₁ always leads to a higher density of particles in region 1, irrespective of whether a₁,a₂ ≫ 1 or not (see figure 3(a)). A physical explanation may be that first, a higher a leads to a lower accumulation of particles at the wall (equation (4)); and second, the exponent (i.e. the eigenvalues) scales with a (equation (25) and so a higher a will lead to more rapid decay of density from the wall (figure 3(c)). These two complimentary effects turn out to be dominant in determining the gradient of the particle density, and thus lead to the observation that a higher a always leads to a lower density of particle in that region. Of course, the validity of this argument is only confirmed by previous analytical investigation.

In summary, I have studied the confinement-induced aggregation phenomenon in a system of self-propelled particles inside a channel. The model system under consideration is drastically simplified in order to focus on the effects of confinement on aggregation. In particular, the particles are assumed to be non-interacting and the walls are assumed to be perfectly rigid and smooth. Even in this simplified scenario, solving the rotationally continuous model is mathematically difficult and thus forbids an analytical understanding of the system. Therefore, I introduced a further-reduced model where the particles can only take six different directions of travel. This enabled me to solve the model analytically and to prove that when the rotational fluctuation a is uniform in the system, the ratio of the number of particles inside the channel versus that accumulated at the walls increases linearly with a in both the a ≪ 1 and a ≫ 1 regimes. I also demonstrated that if one side of the channel has a larger rotational fluctuation, this inevitably leads a lower density of particles on that side. Given the fundamental nature of the model and the general nature of the observations, the effects of confinement on active systems discussed here should continue to be present in more elaborate models and in real experimental systems. In particular, it would be interesting to study to what extend confinement-induced gradient formation in the bulk of an active system can lead to the emergence of pattern formation in the whole system.