Boundary induced convection in a collection of polar self-propelled particles

TL;DR

This study investigates how boundary conditions influence convection patterns in a confined system of polar self-propelled particles, revealing the formation of rolls akin to Rayleigh-Benard convection.

Contribution

It introduces a hydrodynamic model for polar self-propelled particles with specific boundary conditions, showing how density inhomogeneities and convection rolls emerge.

Findings

Density inhomogeneities increase with self-propulsion speed.

Convection rolls form similar to Rayleigh-Benard convection.

Density remains homogeneous at zero self-propulsion speed.

Abstract

We study a collection of polar self-propelled particles confined to a long two-dimensional channel. We write the coupled hydrodynamic equations of motion for density and polarisation order parameter. At two confined boundaries, density is fixed to the mean and orientation is anti-parallel with fixed magnitude of polarisation. Such boundary conditions make our system similar to a sheared suspension of self-propelled particles, which has many practical applications. Antiparallel alignment at the two confined boundaries and alignment inside the channel create rolls of orientation along the long axis of the channel. For zero self-propulsion speed, density and orientation fields are decoupled and density remains homogeneous inside the channel. For finite self-propelled speed, density inhomogeneities develop and these rolls move along the long axis of the channel. Density inhomogeneity increases sharply with increasing the self propulsion speed and then reaches a maximum and again decreases for very large speeds. Formation of rolls is very similar to the classic problem of Rayleigh-Benard convection in fluid dynamics.

Equations

$$\frac{\partial \rho }{\partial t}=-v_0\nabla \cdot (\rho \mathbf{P})+D_{\rho }{\nabla }^2\rho$$

$$\frac{\partial \mathbf{P}}{\partial t}$$

$$r^{\prime }=r/l;$$

$$t^{\prime }=t{D}_R.$$

$$\frac{\partial \rho }{\partial t}=-R_A\nabla \cdot (\rho \mathbf{P})+{\bar{D}}_{\rho }{\nabla }^2\rho$$

$$\frac{\partial \mathbf{P}}{\partial t}=({\alpha }_0-{\alpha }_1|P{|}^2)\mathbf{P}-\frac{R_A}{2{\rho }_0}\nabla \rho +{\bar{D}}_P{\nabla }^2\mathbf{P}.$$

$$R_A=\frac{v_0}{l{D}_R}$$

$${\bar{D}}_{\rho ,P}=\frac{D_{\rho ,P}}{l^2{D}_R}.$$

$$\frac{\partial \rho }{\partial t}=-R_A\nabla \cdot (\rho \mathbf{P})+{\bar{D}}_{\rho }{\nabla }^2\rho +\nabla \cdot {\mathbf{f}}_{\rho }(\mathbf{r},t)$$

$$\frac{\partial \mathbf{P}}{\partial t}=-(-{\alpha }_0+{\alpha }_1|P{|}^2)\mathbf{P}-\frac{R_A}{2{\rho }_0}\nabla \rho +D_P{\nabla }^2\mathbf{P}+{\mathbf{f}}_{\mathbf{P}}(\mathbf{r},t)$$

$$<f_{\rho }^i(\mathbf{r},t)f_{\rho }^j({\mathbf{r}}^{\prime },t^{\prime })>=2{\Delta }_{\rho }{\delta }_{ij}\delta (\mathbf{r}-{\mathbf{r}}^{\prime })\delta (t-t^{\prime })$$

$$<f_{\mathbf{P}}^i(\mathbf{r},t)f_{\mathbf{P}}^j({\mathbf{r}}^{\prime },t^{\prime })>=2{\Delta }_P{\delta }_{ij}\delta (\mathbf{r}-{\mathbf{r}}^{\prime })\delta (t-t^{\prime })$$

$$\Delta (t)$$

$$\displaystyle+|r^{n_{e}}_{min}(i,t_{0})-r^{n_{e}}_{min}(i,t+t_{0})|^{2}\big{]}$$

$$D(t)=\frac{1}{2t}\Delta (t)$$

Full text

I Introduction

Collective behaviour of active particles are extensively studied in sriramrev1 ; sriramrev2 ; vicsekrev ; sriramrev3 . Large collections of living organisms are known to exhibit highly coherent collective motion vicsek1995 ; vicsek2 ; tonertu ; traffic ; shradhapre This behavior, often referred to as “flocking” spans an enormous range of length scales and is seen in diverse systems animalgroup ; helbing ; feder ; kuusela31 ; hubbard ; rauch ; benjacob ; harada ; nedelec ; schaller ; vnarayan . These systems are regoursly studied in bulk either (i) using hydrodynamic equations of motion for slow variables (ii) or microscopic rule based models viz.: Vicsek’s model vicsek1995 . But most biological systems are confined to thin geometry examples . Confinement and boundary plays an important role in variety of problems in biology examples , sheared systems sheared and other places like in fluid dynamics. One classic example include Rayleigh-Benard (RB) convection in fluid rb . In these confined systems, the effect of boundaries are very important.

Boundary can play very important role in a collection of self-propelled particles. It can induce many interesting phenomena like, in many cases, boundary can induce spontaneous flow inside the channel faradaydiscussion . We write the phenomenological equations of motion for local density and polarisation order parameter for the collection of polar self-propelled particles Eqs. 1 and 2. Self-propelled speed (SPS) of the particle introduces a non-equilibrium coupling between density and polarisation. For zero SPS both density and polarisation are decoupled. We solve these equations in the confined geometry shown in Fig. 2. At the two boundaries of the channel orientation of rods are antiparallel, which produces a gradient along the confinement direction. Diffusion tries to make them parallel. Hence the competition between above two create rolls of orientation along the long-axis of the channel. For zero SPS these rolls are static and density inside the channel is homogeneous. For non-zero SPS both density and polarisation are coupled and such coupling produces moving rolls.

In Fig. 1, we show the (left) vector plot of orientation and (right) density of particles inside the channel for SPS $v_{0}=1.5$ or activity $R_{A}=0.67$‘ at different times. We find inhomogeneous moving pattern of orientation and density along the long axis of the channel Fig. 1(top to bottom). Arrow indicate the direction of motion.

In rest of the article, section II discusses the model in detail. Here we also write the hydrodynamic equations of motion for density and polarisation. Section III discusses the numerical details for solving these equations. We discuss our results in section IV and finally conclude with discussion and future aspect of this study in section V.

II Model

We consider a collection of self-propelled particles of length $l$ confined to a two-dimensional channel whose thickness $d$ is very small compare to its long axis $L$. We fix the length of the channel $L$ and vary the width of the channel $d<<L$. Orientation at the lower boundary is parallel to horizontal axis and at the upper boundary it is antiparallel and magnitude of polarisation fixed at two boundaries. We also maintain mean density at two confined boundaries to avoid accumulation of particles at boundaries. Periodic boundary condition is used for both density and polarisation along the long axis of the channel. Geometry of confined channel and orientation of particles at the two boundaries is shown in Fig. 2.

II.1 Hydrodynamic equations of motion

Dynamics of the system is described by the equations of motion for hydrodynamic variables for the collection of polar self-propelled particles. We write the phenomenological coupled hydrodynamic equations of motion for density $\rho$, because total number of particles are conserved and polarisation $P$, which is an orientation order parameter, is a broken symmetry variable in the ordered state. We write the minimum order terms allowed by symmetry. Two equations are

[TABLE]

and

[TABLE]

Density equation 1 is a continuity equation because total number of particles are conserved. Right hand side of Eq. 1 can be written as a divergence of a current ${\bf J}$ which has two parts, ${\bf J}={\bf J}_{a}+{\bf J}_{p}$, where ${\bf J}_{p}\propto\nabla\rho$ is proportional to gradient in density we call it “passive current” and ${\bf J}_{a}\propto v_{0}\rho{\bf P}$ is proportional to polarisation vector ${\bf P}$ and self-propelled speed $v_{0}$ and we call it as “active current”. For zero SPS $v_{0}=0$ or polarisation ${\bf P}=0$ active current is zero.

The order parameter equation 2 contains (i) mean field order disorder terms $\alpha_{0}$ and $\alpha_{1}$, (ii) pressure term present because of density fluctuation and (iii) diffusion in polarisation $D_{P}$. $\alpha_{0}$ and $\alpha_{1}$ are positive and determine the mean field value of polarisation ${\bf P}$ in the bulk $|P|=\sqrt{\frac{\alpha_{0}}{\alpha_{1}}}$. We choose $\alpha_{0}=\alpha_{1}=1.0$, such choice of $\alpha_{0}$ and $\alpha_{1}$ prefers homogeneous polarised steady state $|{\bf P}|=1.0$ in the bulk. $\nabla\rho$ is the pressure term and proportional to the self-propelled speed $v_{0}$ of the particle. $D_{P}$ term is written in the limit of equal elastic constant approximation for splay and bend deformations in two dimensions and $D_{R}$ is the rotational diffusion.

We rescale all lengths by particle length $l$ (which we choose $1$) and time by rotational diffusion time $D_{R}^{-1}$.

[TABLE]

and write the dimensionless equations of motion for density

[TABLE]

and polarisation order parameter

[TABLE]

where

[TABLE]

is the dimensionless activity. It is a ratio between the self-propelled speed $v_{0}$ and the rotational diffusion $D_{R}$. Hence we can increase activity either by increasing speed $v_{0}$ or having small $D_{R}$. We can also define dimensionless diffusions.

[TABLE]

which is again ratio between bulk diffusion and rotational diffusion.

III Numerical Study

We numerically solve the coupled hydrodynamic equations of motion for density and polarisation order parameter. We go beyond mean field and add Gaussian random white noise to the density

[TABLE]

and the order parameter

[TABLE]

random forces are chosen to have zero mean and correlations

[TABLE]

and

[TABLE]

where $\Delta_{\rho}$ and $\Delta_{P}$ are dimensionless noise strengths. Numerical study is done for fix noise strength $\Delta_{\rho}=\Delta_{P}=0.05$. We fix $D_{R}=0.1$, mean density $\rho_{0}=0.1$, diffusivities $D_{\rho}=D_{P}=1.0$. Hence activity $R_{A}$ is varied by changing the SPS $v_{0}$. Typical diffusive length scale in the system $\delta=\sqrt{\frac{D_{P}}{D_{R}}}$ and for the above specific choice of parameters $\delta=\sqrt{10}\simeq 3.5$. Hence, we choose the lower limit on $d>\delta$ such that $d$ is not much bigger than $\delta$ and upper limit such that effect of confinement is important. Hence $d$ is varied from $5$ to $15$. We fixed the length of the channel and SPS is changed from from zero to large values.

We solve these PDE’s 8 and 9 using Euler method for numerical differentiation on a grid $\Delta x=1.0$ and $\Delta t=0.1$ (we have checked that numerical scheme is convergent and stable for the above choice). Inside the channel we start from initially random order parameter and homogeneous density $\rho=\rho_{0}=0.1\pm 0.01$, At the two boundaries the magnitude of polarisation is fixed to $P_{0}=1.0$ and antiparallel orientation such that $P_{x}(z=1)=1.0$, $P_{z}(z=1)=0.0$ and $P_{x}(z=d)=-1.0$, $P_{z}(z=d)=0.0$ and density is maintained to mean density $\rho(z=1)=\rho(z=d)=\rho_{0}$. Periodic boundary condition is used along the long axis of the channel.

IV Results

We divide our results in two parts. First we discuss what happen when activity is turned off, and then we discuss the effect of activity.

IV.1 Zero activity $R_{A}=0$ or SPS $v_{0}=0.0$

For zero activity $R_{A}=0.0$ density and order parameter are decoupled Eqs. 8 and 9. Steady state solution for density is homogeneous density $\rho=\rho_{0}$. Orientation order parameter shows formation of rolls along the long axis of the channel. These rolls are formed because of competition between antiparallel boundary condition which produces a gradient in orientation and diffusion term which tries to make them uniform. Mechanism is similar to Rayleigh-Benard convection in fluid dynamics, where competition between temperature gradient and gravity produces convective rolls. In Fig. 3 (a) we plot the vector plot of orientation field for zero self-propelled speed $v_{0}=0.0$. We find formation of periodic pattern or rolls along the long axis of the channel. In Fig. 3(b) and (c) we also plot the $x$ and $z$ component of orientation field along the long axis of the channel which confirms the above periodic structure. In Fig. 3 (d) we plot the density along the long axis of the channel, which remains homogeneous inside the channel.

IV.2 Non-zero activity or SPS $v_{0}\neq 0.0$

When we switch on the activity parameter $R_{A}$, hydrodynamic equations of motion for density and polarisation order parameter Eqs. 8 and 9 are coupled. Non-zero activity introduces an active current ${\bf J}_{a}\propto v_{0}{\bf P}\rho$ which is proportional to the SPS $v_{0}$ and local polarisation ${\bf P}$ (discussed in detail in section II). This active current produces density inhomogeneity and hence enhanced pressure in local polarisation, which is again proportional to SPS $v_{0}$ or activity parameter $R_{A}$ as shown in Eq. 8 and 9.

IV.2.1 Density fluctuation

We first calculate the density inhomogeneity for different width of the channel. Density inhomogeneity increases as we increase the width of the channel. In figure 4 we plot the percentage density fluctuation $\%\Delta\rho$ along the long axis of the channel averaged over transverse direction for five different widths $d=5$, $7$, $9$, $11$, $13$ as a function of activity parameter $R_{A}$. For very small width $d=5$, density fluctuation is small (strong confinement). For width of the channel $d\geq 7$, density fluctuation shows non-monotonic behaviour as a function of activity $R_{A}$. As we change $R_{A}$ from zero, first increases very sharply with a peak at some finite $R_{A}$ and then decreases slowly for larger activity. Peak position shift towards smaller activity as we increase width of the channel. It changes from $R_{A}=1.13$ to $0.4$ as we change the width from $d=7$ to $d=13$. Hence confinement suppresses the large density fluctuation present in general in self-propelled particles in bulk. Suppression of density fluctuation because of confinement is is also found previously in the study of sheared suspension of Self-propelled particles sheared .

Non-monotonic nature of curve gives a finite range of activity, where fluctuations are large. For very small activity coupling of density to background periodic rolls is small and hence small density fluctuation. As we increase activity, active contribution to density current increases: which is proportional to the local order parameter ${\bf P}$ and hence background periodic rolls of polarisation order parameter. Activity plays two role here, (i) it produces an active density current proportional to local polarisation. Since such active current creates density inhomogeneity hence (ii) it creates pressure in local polarisation. Hence for very large activity density current will be large but at the same time pressure will also increase and it will destroy the background periodic rolls. Hence small active current (because ${\bf J}_{a}\propto{\bf P}$), hence small density inhomogeneity.

IV.2.2 Travelling rolls

Even for zero activity as discussed in section IV.1, antiparallel boundary conditions at the two confined boundaries creates rolls of orientation field. For the range of activity when density fluctuation is large, these rolls move from one end to another end of the channel. We call them as travelling rolls: where density and orientation both shows periodic pattern (please see fig. 1 for one such activity $R_{A}=0.67$). In Fig. 5 we plot density for width of the channel $d=13$, along the long axis of the channel averaged over transverse direction. We calculate density for four different times (with equal time difference) and for four different activity strength $R_{A}=0.13,0.53,1.13$ and $1.87$. For all activities density shows inhomogeneous periodic pattern. Density inhomogeneity shows variation for different activity. For very small activity $R_{A}=0.13$, density shows small inhomogeneity and remains static with time. For activity $R_{A}=0.53$, density is periodic as well as inhomogeneous. These periodic pattern move with time. In Fig. 5(b) we draw a horizontal arrow to denote the motion of periodic pattern with time. Larger the arrow faster the pattern move. Direction of arrow shows the direction of motion of periodic pattern. This direction is spontaneously chosen from two equal possible directions in the system. For activity $R_{A}=1.13$ we get the similar result as for $R_{A}=0.53$ but density inhomogeneity is weaker. For activity $R_{A}=1.87$, density is periodic but inhomogeneity is even more weak and very small shift in peak position as a function of time (as shown by small horizontal arrow in Fig. 5(d)). Hence travelling periodic pattern and density inhomogeneity are coupled and larger density inhomogeneity creates faster moving rolls from one end to other end of the channel. For very large activity, density coupling is strong enough such that it destroys the background periodic pattern of orientation and hence no travelling rolls.

IV.2.3 Mean square displacement

We further characterise properties of travelling rolls assuming periodic density pattern as shown in Fig. 5. Cartoon picture of maxima and minima of periodic density profile for one such realization is shown in Fig. 6. We track the position of maxima and minima of periodic density pattern. Position of each maxima and minima we model as position of independent particles moving in one dimension. Hence each maxima and minima represent one particle and we save the position of maxima and minima or particles position with time. We calculate the square displacement of such positions and take average over all maximas and minimas and many initial realisations. Hence mean square displacement is defined as

[TABLE]

where $r^{n_{e}}_{max/min}(i,t)$ is the position of maxima/minima of the $i^{th}$ periodic profile at time $t$ for $n_{e}^{th}$ realisation (as shown in Fig. 6). Averaging is done over all periodic positions $i=1,N_{i}$ and number of realisations $n_{e}=1,N_{e}$. For our calculation we used total number of realisation $N_{e}=20$.

When travelling rolls form as discussed in previous section IV.2.2, then $\Delta(t)$ is proportional to $t^{2}$. In Fig. 7(a) we plot MSD, $\Delta(t)$ for different activity $R_{A}$ for channel of width $d=13$. For zero activity density is homogeneous and no periodic pattern and $\Delta(t)\simeq t$ (diffusive). As we increase activity $R_{A}=0.13$, MSD $\Delta(t)$ is subdiffusion where $\Delta(t)\simeq t^{\alpha}$ and $\alpha<1$. Subdiffusive behaviour shows the arrest of density in the center of the periodic rolls, which acts like a disorder site (with small polarisation). Similar arrest of density in the presence of quenched disorder field is found in recent study of peruaniprl . But in our model disorder site or center of the roll moves for sufficient large activity. For activity $R_{A}\geq 0.67$, MSD shows two regimes with initial Subdiffusive with $\alpha<1$ to later travelling motion with $\alpha\simeq 2$. Hence for $R_{A}\geq 0.67$, density remains arrested for some time and then travelling rolls sets in. As we further increase activity $R_{A}=0.67$ and $1.33$ we find initial Subdiffusive with $\alpha<1$ and the later travelling motion with $\alpha=2$. Time spent in arrested state decreases as we increase activity. For large activity $R_{A}\geq 1.67$, $\Delta(t)\simeq t^{\alpha}$ for very small time and then switches to transient faster dynamics and then saturates to diffusive $\Delta(t)\simeq t$ for large time.

In Fig. 7(b) we also plot the diffusivity defined as

[TABLE]

$D(t)$ remain flat for zero activity, hence diffusion. For small activity $R_{A}=0.13$ it decrease with time and shows arrested subdiffusion. For intermediate activities $R_{A}=0.67$ to $1.33$, $D(t)$ shows initial subdiffusion and later travelling motion with $D(t)\simeq t$. For very large activity $R_{A}\geq 1.67$ initial $D(t)$ decreases with time then faster growth to diffusive regime with constant $D(t)$ for very long time.

We further explain the dynamics of particle inside the channel as we vary the activity $R_{A}$. For zero activity density and orientation are decoupled hence we expect normal diffusive motion for the particle. As we have discussed before antiparallel boundary condition at the two confined boundaries creates periodic ${\it rolls}$ along the long axis of the channel. As we switch on activity density is coupled to periodic orientation field. For very small activity particles are trapped in periodic orientation. center of periodic pattern acts as quenched disorder site. For intermediate activity coupling is strong and large density inhomogeneity as shown in Fig. 4. Hence for initial transient time particles are trapped to periodic pattern of orientation and later active current sets the travelling rolls. Such initial subdiffusion to propagating motion is very common for particle moving in periodic structure. When particle move in a periodic background it spent some of its time in trapped phase and then start propagation. Similar subdiffusion to propagation is found for the particle moving in periodic media shradhajstatphys . For very large activity coupling is very strong and it destroys the background periodic pattern of orientation and hence diffusive behaviour at late time.

V Discussion

In our present work we write the phenomenological hydrodynamics equations of motion for density and local polarisation order parameter Eq. 8 and 9 for the collection of self-propelled particles. We solve these equation in the confined channel of width $d\ll L$ very small compare to the long axis of the channel. We impose an antiparallel boundary conditions at the two confinement boundaries and maintain the density at its mean value and the magnitude of polarisation constant. Such a geometry is important because it mimics the shear. First we solved equations(1,2), for zero activity. Antiparallel boundary conditions impose a gradient of orientation along the confinement direction and diffusion tries to make them parallel. Hence there is a competition between these two terms and we find periodic patterns of the orientation field along the long axis of the channel Fig. 3(a). And, since for zero activity density is not coupled to the orientation field, it remains uniform.

Non-zero activity, turn on a contribution of active current, which is proportional to the local polarisation inside the channel. Such active currents make the inhomogeneous density inside the channel. Density inhomogeneity increases as we increase the width of the channel. For fixed channel width, as we increase activity initially density inhomogeneity increases with activity $R_{A}$ Fig. 4 and then decrease for large activity $R_{A}$. For fixed channel width, for the range of activity when density inhomogeneity is large, density periodic pattern sets in and start moving from one end to other end of the channel. Real space image of moving periodic orientation rolls and density profile for fixed channel width and activity is shown is Fig. 1 for different times.

Travelling periodic rolls we observe here is very similar to Rayleigh-Benard convection in fluids. Orientation plays the role of temperature gradient and diffusion is like gravity which acts opposite to gradient. A competition between these two produces periodic rolls and in the presence of activity these rolls move from one end to other end of the channel.

It would be interesting to study other kinds of boundary condition on the flow properties of active particles inside the channel. For example recent study of ncomm , where boundary induces accumulation of particles. This phenomena where boundary induces spontaneous flow in confined channel can give some insight of transport of active fluid in biology biologicaltransport .

Acknowledgements.

S. Mishra would like to thank DST India for financial support. S. Mishra gratefully acknowledge helpful conversation with Abhik Basu and Argha Banerjee.