Dynamics of sedimenting active Brownian particles

We investigate the stochastic dynamics of one sedimenting active Brownian particle in three dimensions under the influence of gravity and passive fluctuations in the translational and rotational motion. We present an analytical solution of the Fokker--Planck equation for the stochastic process which allows us to describe the dynamics of the active Brownian particle in three dimensions. We address the time evolution of the monopole, the polarization, and the steady-state solution. We also perform Brownian dynamics simulations and study the effect of the activity of the particles on their collective motion. These results qualitatively agree with our model. Finally, we compare our results with experiments [J. Palacci \emph{et al.}, Phys. Rev. Lett. \textbf{105}, 088304 (2010)] and find very good agreement.


Introduction
Active particles convert energy from chemical, biological, or other processes into motion. The study of active particles, and especially their collective motion, has received much attention due to a renewed interest in the physical principles underlying the motion of, e.g., plankton or bacteria, and also on account of technological applications involving both biological and artificial controllable active systems [1][2][3]. Active particles exhibit a fascinating multitude of interesting behaviors from the single particle to collective states [4][5][6], due to their nonequilibrium nature.
Typically, active particles move in an aqueous environment, where, because of their size, viscous forces dominate, and inertial forces are completely negligible. In fact, the consideration of the Navier-Stokes equations identifies that the nature of the dynamics is dictated by the ratio of viscous to inertial forces, known as the Reynolds number R = σvρ/η, where σ is the typical size of the microorganism, v its mean velocity, and ρ, η are the fluid's density and viscosity, respectively. For motile bacteria R≈10 −5 . As noted by Purcell [7], this means that if the propulsion of the active particle were to suddenly disappear, it would only coast for 0.1Å. Thus, the state of motion is only determined by the forces acting at that very moment, and inertia is negligible.
Even in dilute suspensions, where particle-particle interactions can largely be neglected, and the dynamics are a e-mail: jeremy.vachier@ds.mpg.de b e-mail: marco.mazza@ds.mpg.de dominated by the balance of active motion and gravity, interesting results are found [8][9][10][11]. Palacci et al. [12] showed experimentally with active Janus colloids that activity increases the sedimentation length, by increasing the effective diffusivity. More recently, Ginot et al. [13] characterized the equation of state of sedimenting active colloids as a function of the activity.
Theoretical studies of active particles, based on the framework of active Brownian particles [14] and stochastic processes, have mostly focused on two-dimensional systems [15][16][17][18][19]. A complete description in three dimensions (3D) in terms of the Fokker-Planck equation is challenging [20,21] and some recent progress in the theory of one active particle [9,22,23] highlights the fact that many questions are still open, especially in 3D. For example, in dilute suspensions, what is the transient sedimenting dynamics? the emergence of polarization (and possibly higher orders) is intriguing and currently under investigation [24]; what are the appropriate variables to construct an equation of state? In denser suspension, the important role of hydrodynamic interactions makes the situation even more complicated. For what physical conditions is the sedimenting steady state stable? What are the other possible steady states? Can we write an equation of state in this case? What are its relevant dynamical variables? In this work, we address the first question, that is, the transient state.
We aim to analytically characterize the sedimentation of one active Brownian particle in 3D and, by means of Brownian dynamics simulations for many weaklyinteracting particles. First, we analytically describe the sedimentation of one active particle under gravity with two overdamped Langevin equations and the associated Fokker-Planck equation to obtain the particles' density profile in the direction of gravity. The density profile is obtained from the probability density function P (r, e,t|r 0 , e 0 ,t 0 ) of finding an active particle at the position r, with an orientation e at time t, given the initial state (r 0 , e 0 ,t 0 ). Due to the complexity of the problem, finding the general expression of P (r, e,t|r 0 , e 0 ,t 0 ) in 3D is challenging. The method presented in this article allows the coupling between the orientation and the position obtained in 3D, which we then specialize in one direction. Furthermore, in comparison with previous work [9,19,20,22,25] this method has the additional advantage of providing access to the full temporal dynamics, and is not limited to steady-state conditions. We find an approximate solution for the time-dependent density, polarization and the steady-state solution. Secondly, we perform Brownian dynamics simulations to describe the sedimentation of many particles, where fluid-mediated hydrodynamic interactions are approximated via a short-range potential with up-down symmetry.
The remainder of this work is organized as follows. In sect. 2, we introduce the stochastic process and solve the associated Fokker-Planck equation for a single, sedimenting active Brownian particle. In sect. 3, we show the results of Brownian dynamics simulations of dilute suspensions of active particles. Finally, in sect. 4 we discuss our conclusions.
2 Analytical solution for a single active Brownian particle We study analytically the motion of one self-propelled microscopic particle (active particle), considered as a point particle, in 3D under an external force: gravity. An example of this motion is shown in fig. 1. An active particle has the ability to convert energy into motion. We represent the self-propulsion with a constant speed v s acting on the particle. Typically, an active particle moves inside a fluid and due to its microscopic size we cannot neglect the influence of thermal fluctuations caused by the surrounding fluid buffeting the particle. The interactions with the fluid are represented by stochastic terms as for a Brownian particle. Due to gravity, the suspended active particle approaches a stationary state where its position has an increased probability of being close to the confining interface. This phenomenon is called sedimentation. To describe this motion, we derive a Fokker-Planck equation [26][27][28][29][30].
We treat the motion of one active particle, described as a point particle, in 3D under gravity by considering a single active Brownian particle moving with a constant active speed v s along a direction represented by the orientation e, subject to random fluctuations. The motion of the active particle is biased by a drift velocity −v g z in the direction of gravity. Our system is then described by two overdamped Langevin equations The random fluctuations are modeled in terms of the vectors ξ and ξ e , with zero-mean, Gaussian white noise components, and with variance ξ where D t and D e are the translational and rotational diffusivities, respectively, δ ij is the Kronecker delta, and δ(t) the Dirac distribution. From eqs. (1), (2) we can derive a Fokker-Planck equation that accounts for the evolution in time of the oneparticle probability density function, P (r, e,t|za, e 0 ,t 0 ), of finding an active particle under gravity diffusing in 3D, with the initial condition P (r, e,t = t 0 |za, e 0 ,t 0 )= δ(r−za)δ(e−e 0 ). In the following, to lighten the notation we will use P (r, e,t)=P (r, e,t|za, e 0 ,t 0 ).
After some manipulation (see appendix A), we obtain the following Fokker-Planck equation +D t ∇ 2 P (r, e,t)+D e L e P (r, e,t), with ∂φ 2 ] the Laplace-Beltrami operator on the 2-sphere S 2 , and where we expressed e =( s i nθ cos φ, sin θ sin φ, cos θ) T in spherical coordinates. Equation (3) can be written symbolically as a continuity equation which defines the current J . To solve eq. (3), we proceed in the following way: i) we will move to Fourier space; ii) we will use an expansion in terms of eigenfunctions of the Fokker-Planck operator; iii) we will perform a multipole expansion; iv) we will focus on the dependence of the probability on the z-direction, along which gravity applies; and v) we will perform the inverse-Fourier transform. The Fourier transform of eq. (3) then reads ∂ ∂t P (k, e,t)=iv s e · k P (k, e,t) − iv g k z P (k, e,t) −D t k 2 P (k, e,t)+D e L e P (k, e,t).
Let us for the moment consider the simple case v s =0a n dv g = 0 which corresponds to a simple Brownian particle. Because the operator is Hermitian, its eigenfunctions e −Dtk 2 t e −λnDet Y m n (θ, φ) form an orthonormal basis of the space of our solutions, where λ n = n(n + 1), (n−m)! (n+m)! P m n (cos(θ))e imφ are the spherical harmonics including the Condon-Shortley phase factor, and n, m ∈ N, −n ≤ m ≤ n. Additionally, because O g = iv g k z is simply a multiplicative scalar, an eigenfunc-

Taking into account the initial condition
where the coefficients P m n (k,t) are determined by imposing that the expression in eq. (6) satisfy eq. (5) (see appendix B). Physically, the infinite sums on the right-hand side of eq. (6) represent the increasingly faster decay with time of higher-order spherical harmonics [21].
Because of the rotational dynamics in our problem, it is convenient to explicitly highlight the underlying physical symmetries by expanding the full probability P (r, e,t) in terms of spherical tensors, that is, the irreducible representations of the rotation operator. Each spherical tensor transforms like the eigenfunctions of the angular momentum of corresponding rank n =0 , 1, 2,..., where the first three tensors represent the density ρ (monopole, n = 0), the polarization D (dipole, n = 1), and the nematic tensor Q (quadrupole, n = 2), respectively. The probability can then be expanded as P (r, e,t)=ρ(r,t)+D(r,t) · e + e · Q(r,t) · e + .... (7) In the large time limit, the monopole term ρ(r,t) will dominate the sedimentation process (while higher-order terms in eq. (7) are relevant for observables with shorter characteristic time scales). Its Fourier transform can be found by truncating the sum in eq. (6) at the n =0term, that isρ After some computations (see appendix C), we can work out the equation governing the dynamics of P 0 which is the telegrapher's equation [31], and accounts for processes with a finite speed of propagation. A solution of eq. (9) reads P 0 a n d G(k z )a rbitrary functions of the wavevector in the z-direction k z . The expression for the monopole is found from the inverse Fourier transform, and reads in the exponential makes it difficult to perform the inverse Fourier transform. However, because we are interested in the long-wavelength limit of the sedimentation profile, it is natural to consider a Taylor expansion of w(k z ) around k z =0 where F and G are defined as follows and similarly for G. Elementary integration yields where we have defined the effective diffusivities D ± eff ≡ D t ± repeatedly observed in experimental [12] and theoretical works [23,32]. By imposing that mass is conserved during the sedimentation process we can determine the functions F =exp(2D e t)and G =1.
In order to describe the sedimentation phenomena, we need to impose a reflective boundary condition at the confining wall located at z = 0. By integrating eq. (3) over the orientation, we find The associated continuity equation reads and by imposing no flux J z = 0 at the wall, the boundary condition, of the Robin type, reads [33] =0. (11) Moreover, in large time limit, the monopole ρ(z, t) dominate the sedimentation process and therefore the dipole We rewrite eq. (10) as and In order to impose no net flux across the reflective wall, we use the method of images, but as known in the theory of partial differential equation (Sommerfeld) the appropriate image system consists of replacing the wall at z = 0 with a mirror source placed at z = −a (in addition to the real source at z = a) and a continuous sequence of images which take place at all points ξ<−a [34][35][36][37][38]. We can rewrite the probability density as By applying eq. (13) to our system, ρ r (z, t) reads where the coefficients A 1 , A 2 , k 1 (ξ), and k 2 (ξ) are also found via the Robin boundary condition (see appendix D) and the solution yields The steady-state regime is given by taking the limit t →∞ of eq. (15) and reads (16) We note that the correct boundary condition is given by eq. (11), but in the large time limit the dipole term vanishes; additionally, by matching our parameters with the experimental values in [12], for which D − eff is positive, it follows that D t >v 2 s /6D e and D 2 t ≫ v 4 s /36D 2 e , and that the Robin condition is satisfied in the large time limit.
In the following, we take the active particle's diameter σ, mass m and its translational diffusion coefficient D t as the units of length, mass and diffusivity. Thus, we can measure rotational diffusivity in terms of D e = D t /σ 2 .A dimensionless measure of the relative strength of the selfpropulsion to the diffusive behavior, that is, the relative persistence of the active motion, is given by the Péclet number P = v s σ/D t . Figure 2 shows the evolution of the density profile ρ r (z, t) found from our solution to the sedimentation process in eq. (15). The initial position of the active particle is chosen at z/σ = 40. The corresponding initial density ρ r (z, t = 0) is a Dirac delta distribution. As time progresses, we observe the shifting and flattening of the density profile. The steady-state regime, given by eq. (16) is characterized by an exponential decay, which matches the sedimentation profile. We match our parameters with the experimental values given in [12], where the Péclet number 0.5 P 5, and we find a near-quantitative We observe a good match with the experimental results in [12]. The model parameters are Deσ 2 /Dt =1.8, and vs/vg ∈ [1.2, 3.5]. agreement with the experiments. As predicted in theoretical works [20,22], ρ r (z, t) decays exponentially away from the confining surface. Upon increasing the self-propulsion v s , and therefore the effective diffusivity D eff , the density profile tends to spread away from the wall as observed in the experiment [12]. This behavior is shown in fig. 3.
The sedimentation length δ eff is the characteristic length scale of the decay of ρ r (z, t) with z.Itw asfoundto depend strongly on the activity of the self-propelling particle [12,20]. In general, we find a linear relationship governing the growth of δ eff with D eff /v g , δ eff = c 0 +D eff /v g .The constant c 0 ≡ c 0 (v g ), and can be chosen to be zero, which is the value consistent with the experiments in [12]. The relationship between δ eff and D eff provides a connection between the microscopic behavior of the active particle and the long-time emergent dynamics [12].
Moreover, the precise nature of the density profile in proximity of the confining surface will be affected by a number of effects such as: electrostatics, and hydrodynamic interactions of the active particles with the walls. For example, in a recent work [39], the authors show that boundaries can steer Janus colloids, which, as a result, move above the boundary at a fixed distance. These effects are not taken into account here.
Additional information about the active sedimentation process can be gained by considering the next term in the expansion eq. (7), i.e. the polarization. The probability density function becomes P (z, cos(θ),t) ≃ ρ(z, t)+D(z, t)cos(θ).
We can express the polarization D by means of the Legendre polynomials. Again, we are only interested in the z-direction. In Fourier space we find After some computations and applications of the boundary conditions (at z =0,J =0andθ = π), the probability density function reads We refer the readers to appendix E for further details, and for the definition of the constants B 1 , B 2 , C, α,a n df . After imposing the condition P (z, cos(θ),t) ≥ 0, we show in fig. 4 P (z, cos(θ),t), at large t and for three values of the orientation θ ∈{0,π/2,π}. As predicted in [9], we observe an accumulation of active particles with a net polarization at the bottom wall, moving againts the wall (θ = π). The qualitative picture is in agreement with [9].

Simulations of the collective motion
We next investigate a many-particle system composed of active Brownian particles under the effect of gravity. An analytical approach is a formidable task; thus we turn to numerical simulations. Even in the simpler case of passive Brownian particles, sedimentation is a complex process on account of velocity correlations and hydrodynamic interactions [40][41][42]. Because we are interested in the impact of active motion on the sedimentation, we can reduce the nonlinearities associated to hydrodynamic interactions by working in the dilute limit, similarly to [9,12,20,22]. We do however consider weak hydrodynamic interactions via a short-ranged effective potential (see below). We perform Brownian dynamics simulations in 3D.
We consider again a typical Reynolds number R≪ 1. An example of such particles in a biological setting is the microalga Chlamydomonas reinhardtii, which has a typical length σ =1 0μm and self-propulsion speed v s = 60 μms −1 . In colloidal physics active Janus particles have a typical linear size σ =1μm and self-propulsion speed which can vary as a function of the chemical gradient. As a reference, we can take the results found in experiment [12], where the self-propulsion v s =(0.3-4) μms −1 . We describe the system by using two first-order stochastic differential equations, for N active particles where i =1 ,...,N, e i = 1 (implemented by means of a Lagrangian multiplier), v g is the limiting velocity of a particle in the fluid under gravitational acceleration. In eq. (19), φ WCA =4 ǫ[(σ/r ij ) 12 − (σ/r ij ) 6 ]+ǫ is the Weeks-Chandler-Anderson potential [43], r ij = |r i − r j |, representing a hard-core repulsion between active particles, where σ is the linear size of the active particles, and ǫ is the energy scale of the repulsive interaction. In eq. (20), U = i =j cos 2 (θ ij ) is the Lebwohl-Lasher potential, which we use to model to first approximation the updown symmetric interaction due to hydrodynamics that tends to align neighboring particles (see, e.g., [44]). We integrate eqs. (19), (20) with a discretization scheme based on the Euler-Maruyama algorithm in which we take into account the issue of multiplicative noise. We solve eq. (19), (20) in a domain of volume V = L 3 , with L =5 0 σ, with a reflective wall at the bottom at z =0 , and gravity pointing in negative z-direction. The filling fraction of our system is φ = N π 6 σ 3 /V =1 0 −3 . Our results shown below are averaged over 10 4 independent simulations. Figure 5 shows the dependence of the density profile ρ r (z, t) on the position z at different times. At t = 0, the active particles are randomly placed on a plane located at z/σ = 40. After some time, all the active particles sediment on the bottom wall. We observe a qualitative agreement of our simulation with our theory. We conclude that, as long as hydrodynamic interactions are weak, or the system is diluted enough, the theory derived for a single active particle is also applicable to a many-particle system.
To measure the importance of the active motion with respect to the diffusion, we vary the Péclet number P by changing the self-propulsion v s . Figure 6 shows the results of our simulations for the density profile in the long time behavior, and at different values of self-propulsion v s . When the activity is lower than the sedimentation velocity, v s /v g =0.2, we observe a clear sedimentation profile.
If v s /v g =1.0, we observe a weaker sedimentation profile and a peak appears close to the upper part of the sim- ulated domain. This peak is due to balance of the weak, effective hydrodynamic interactions introduced in our simulations with the self-propulsion and the effect of gravity. This is a consequence of the emerging polar order in sedimenting active particles as discussed in sect. 2.
Finally, as soon as v s >v g , we observe an accumulation of particles both on the bottom and top end, which differs from the classical sedimentation profile. We clearly highlight the importance of the activity of the particles, which allows them to move against an external force, in our case the gravity. From a biological point of view, this capacity play an important role, e.g. the algae need light to survive and they need to move against gravity to reach the surface.

Conclusion
The dynamics of sedimenting active particles prove to be an interesting arena where different nonequilibrium effects are at play, providing a testbed for our understanding of far-from-equilibrium phenomena. Including a selfpropulsion to the motion of sedimenting colloidal particles ushers in a wealth of intriguing effects unimaginable from the classical results of Perrin [45]. This is in fact reflected in the considerable interested elicited by this problem [9,12,13,19,20,22].
We studied the sedimentation process of active Brownian particles in three dimensions. Firstly, we developed an analytical method describing the sedimentation profile of one active particle. We solved analytically the Fokker-Planck equation for an active particle in the presence of gravity and a confining wall at the bottom. We addressed the time evolution of the monopole, and found a solution which matches the late-time density profile in [12,13]. Furthermore, we calculated the dipole, and fund the emergence of polar order at the bottom wall, with an accumu-lation of particles moving against the wall, and a depletion of particles moving away from it.
Imposing the no-flux condition at the confining bottom wall produces the steady-state solution. This solution is consistent with a number of previous results (most recently [19], for example).
We recovered the following experimental results found in [12,13]: i) the exponentially decaying density profile for the long-time regime and the steady state ii) the increasing sedimentation length upon increase of the effective diffusivity. Importantly, our method retains the temporal dynamics of the sedimentation process, and therefore in addition to the steady state we also have access to the intermediate states. Our method also allows us to keep the coupling between the rotational and the positional degrees of freedom. In order to characterize more realistic conditions for the sedimentation process, we also considered many particles with weak, effective hydrodynamic interactions and we carried out Brownian dynamics simulations. We are able to measure the importance of the active motion in the long time behavior by varying the Péclet number P. We recovered the density profile, shown in fig. 6, found experimentally in [13] and in numerical simulations [20] of active bottom-heavy particles. However, our model give access to finer details of the sedimentation process of active particles as a function of the activity. Furthermore, the sedimentation profile predicted by our analytical method for one active particle (see fig. 3) matches our simulations for many particles (see fig. 5).

Appendix B. Eigenfunction expansion
Inserting eq. (6) into eq. (5), and then multiplying by the complex conjugate of the spherical harmonics Y * m ′ n ′ and integrating over the solid angle dΩ =sinθ dθ dφ, we find To proceed we define the following integrals: Appendix C. Telegrapher's equation Here we provide details of the computation of the telegrapher's eq. (9). We start by considering the equations for the two coefficients P 0 0 and P 0 combining these two equations yields 3) and after replacing the last term on the right-hand side with eq. (C.2) we find Finally, neglecting the higher order yields which is the telegrapher's eq. (9).

Appendix D. Monopole reflective boundary
Let us start from the solution of a diffusion process It is convenient to introduce a change in the independent variable [38] ρ(z, t)=Ue In order to take into of the reflecting barrier, ρ(r, t)b ecomes [35][36][37][38] ρ r (z, t)=ρ(z, t|a)+Aρ(z, t|−a)+ −a −∞ k(ξ)ρ(z, t|ξ)dξ, (D.1) which tells us that an isolated point (image) z = −a is not sufficient, but we need a continuous sequence of images which take place at all points ξ<−a. The Robin boundary condition for our system reads =0. then By applying the Robin boundary condition (D.1) reads By setting the terms of different time dependence individually equal to zero, the coefficients A and k(ξ) read By replacing the function and coefficient inside the equation, the solution for U r is given by and by rewriting the integral term 2D t e −(z+ξ) 2 /4Dtt dξ and by doing a change of variable [38] for the limits of the integral U r reads Finally, the complete solution reads Simple computations yield where f = 4 4 15 − 2 3 . By applying the inverse Fourier transform D(z, t)= 3 8π 2 C 1 √ 2tD e e −(a−z+vgt+fitDevs) 2 /4tDt 2tD e e −(a−z+vgt−fitDevs) 2 /4tDt .
We only consider the two first terms of the probability density function P (z, cos(θ),t) ≃ ρ(z, t)+D(z, t)cos(θ) and by plugging in the polarization, and focusing only on the real part of the exponentials, the probability density function reads In order to find the coefficients C 1 and C 2 , we apply the Robin boundary conditions at z =0 ,J z = 0, to eq. (4). Moreover, F = e 2Det and G = 1 to keep the mass constant over time. After some computations we find and For the sake of simplicity we set C 2 =1.

Author contributions statement
MGM conceived the project. JV implemented the theory and performed simulations. MGM and JV interpreted the data and wrote the paper. All authors discussed the results and commented on the manuscript.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.