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

I study the confinement-induced aggregation phenomenon in a minimal model of self-propelled particles inside a channel. Starting from first principles, I derive a set of equations that govern the density profile of such a system at the steady-state, and calculate analytically how the aggregation at the walls varies with the physical parameters of the system. I also investigate how the gradient of the particle density varies if the inside of the channel is partitioned into two regions within which the active particles exhibit distinct levels of fluctuations in their directions of travel.


Introduction
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 [1]. Many different physical mechanisms can lead to the aggregation of active systems. For example, these mechanisms can be directional alignment interactions [2,3,4,5,6], density-dependent motility [7,8], and aggregation at the wall when the system is under confinement [9,10,11,12]. 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 [13].
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 (e.g., see [14] 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 [10]. On the other hand, if pairwise interactions between the particles become negligible, such a confinementinduced 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 Active particles are confined in a channel that is infinite in the y-direction but bounded in the x-direction. The walls are assumed to be perfectly rigid in that the particles cannot penetrate them but the particles are free to rotate at and slide along the walls. The channel inside is allowed to have two distinct rotational diffusion constants: a 1 for x ≥ 0 and a 2 for x < 0. (b) The simplified model also considered here restrict the active particles to have only six distinct directions of travel. 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 [15]. To make further analytical process, I will introduce a simplified model in which the particles' directions of travel are discretised.
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 ( Fig. 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 [16].
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.

Model
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 where 1 ≤ i ≤ N, R ≡ (r 1 , . . . , r N ), Θ ≡ (θ 1 , . . . , θ N ), v(θ) ≡ (cos θ, sin θ), and (Eq. 2) is an Itô stochastic differential equation [17] such that The specific geometry considered is depicted in Fig. 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 [18,19]. As explained in [20], this simplification is valid if D ≪ u 2 a −1 , 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.

Intuitive arguments
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, 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 ef f = u 2 /a (see, e.g., Ch. 6 in [21]). On the other hand, if a is small, the motion of the particle inside the channel will be predominantly ballistic. Therefore, Remembering that the units are set such that d = u = 1, τ c becomes Together with (Eq. 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 If the channel is now partitioned into two regions system with distinct rotational diffusion constants a 1 and a 2 ( Fig. 1(a)), an interesting question that is of experimental interests [16] is which side of the region will be of higher density. Focusing on the a 1 , a 2 ≫ 1 regime, one may expect that the particles would behave diffusively inside the channel, with effective diffusion constant D 1 ∼ 1/a 1 and D 2 ∼ 1/a 2 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 localise particles [7,8]. 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 localise certain diffusing proteins [22,23]. If this expectation extends to our active system, then one prediction would be that for a 1 , a 2 ≫ 1, a 1 < a 2 implies n (1) c < n (2) c where n (i) c denotes the number of particles in region i.
Summarising this section, the following three predictions seem to follow from the intuitive arguments presented above: (i) In the symmetric case (a 1 = a 2 = a), if a ≪ 1, then n c /n w ∼ a (?) (ii) In the symmetric case, if a ≫ 1, then n c /n w ∼ a 2 (?) (iii) In the asymmetric case, if a 1 < a 2 and a 1 , a 2 ≫ 1, then n (1) c < n (2) c (?) With the mathematical tools developed in the next section, I will show that in fact only the first "prediction" is true.

Mathematical formulation
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 [24]: 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(r 1 , θ 2 ) ≡ N dr 2 · · · r N dθ 2 · · · θ N f (R, Θ), which is governed by the following dynamical equation [25,26]: Focusing on the steady-state solution and assuming the channel is infinite along the y-direction, the equation simplifies to Employing the separation of variables method by assuming that ρ(x, θ) = A(θ)B(x), I arrive at the two equations where q corresponds to a set of constants to be determined. Eq. (12) can be easily solved: while Eq. (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 [15].
Since the characteristic numbers for these types of Mathieu functions are identical for ±q, the expansion for p is a n e q 2n ax ce 2n (θ/2, 2q 2n ) + b n e −q 2n ax ce 2n (θ/2, −2q 2n ) . As usual, the expansion parameters are obtained by the boundary conditions, which I will determine now.
Similarly, the equation governing the steady-state profile at the wall on the right R(θ) is 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 Fig.  2(a)). To account for these influxes of particles from the walls, Eq. (10) is modified as follows: In summary, the steady-state density profile of the active particles is obtained by solving the coupled differential equations Eqs (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 (Eq. (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 [15,27]. 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.

Simplified 6-direction model
In this simplified model, the angular components are partitioned into six directions as shown in Fig. 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 i-th 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 1 = p 3 , p 6 = p 4 , and also that the functions p only depend on x but not y. As a result, the dynamical equations for the system are ∂ t p 6 = a(p 1 + p 5 − 2p 6 ) + b∂ x p 6 (22) where b = √ 3/2 since b = cos π/3. In matrix form, dp/dx = Mp where p = (p 1 , p 2 , p 5 , p 6 ) T and The four eigenvalues associated to the matrix M are {0, 0, −aγ/b, aγ/b} where γ = √ 3 + 4b + 4b 2 , and the corresponding eigenvectors are the column vectors in the following matrix: Furthermore, the steady-state solution is of the form: where σ k ∈ {0, 0, −aγ/b, aγ/b} are the corresponding eigenvalues of M. The above expression is to be compared to Eq. (14) in the continuum model.
To determine the coefficients c k in Eq. (25), I employ the discrete version of the boundary conditions as discussed in the previous section (Eq. (15)), which for the wall on the left gives, and the flux of particles out of the wall on the left is (see Eq. (16)) where the second equality follows form Eq. (28). Similar expressions can be readily found for the R i . In particular, the condition in Eq. (30) and the corresponding one for p R i enable me to determine the coefficients in the expansion in Eq. (25) completely.

Symmetric case (a 1 = a 2 = a)
In the symmetric case, the interesting quantity to investigate is the ratio of the number of particles inside the channel n c vs. the number of particles at the wall n w . The analytical expression for n c /n w can be obtained straight-forwardly with the help of Mathematica [28]. Although the analytical expression is complicated and difficult to decipher, it does provide the following asymptotic expressions with respect to a: 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 (Fig. 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.

Asymmetric case (a 1 = a 2 )
In the two-region case as depicted in Fig. 1(a), the expansion coefficients in Eq. (25) are obtained by employing the boundary conditions at both walls (Eq. (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 (Fig. 3(c)). Fig. 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 ( Fig. 3(a)) or includes the particles at the walls ( Fig. 3(b)). In other words, if a 1 < a 2 , then n (1) c > n (2) c . Hence, the expectation (iii) discussed in Sect. 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 [16].

Discussion
As we have seen, among the three expectations arisen from intuitive arguments in Sect. 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 Sect. 3, a higher a 2 in relation to a 1 always leads to a higher density of particles in region 1, irrespective of whether a 1 , a 2 ≫ 1 or not (see Fig. 3(a)). A physical explanation may be that first, a higher a leads to a lower accumulation of particles at the wall (Eq. 4); and second, the exponent (i.e., the eigenvalues) scales with a (Eq. (25) and so a higher a will lead to more rapid decay of density from the wall (Fig. 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 vs. 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.