Numerical Simulations of Active Brownian Particles

Active particles are microscopic particles capable to propel themselves, such as bacteria, chemotactic cells, organelles within the cells, and self-phoretic colloids. For their very nature, they belong to the class of non-equilibrium phenomena. They play a significant role as a model system in our understanding of the dynamics of out-of-equilibrium systems, especially those involving living organisms at the microscale. Furthermore, the interest in these active particles is also fueled by their potential for applications as autonomous agents for health care, sustainability, and security applications.


Introduction
Active particles differ from their passive counterparts for their ability to propel themselves. In Fig. 7.1 examples of active particles are given, classified with respect to their size and their propulsion speed. Living microorganisms propel themselves for different purposes such as finding food, escaping from predators or other dangers, and patrolling a territory [1]. Inspired by these microorganisms, researchers have recently developed several artificial particles capable of self-propelled motion activated by localised light, concentration, temperature gradients [2]. Despite the variety of possible self-propulsion mechanisms, we can identify some key features to describe the motion of a self-propelling micro-or nanosized particle: (1) directionality over a characteristic time interval, (2) orientational noise, and (3) absence of inertia. We note that in the case of living organisms the self-propulsion mechanism often implies a deformation in their shape; however, for simplicity, this aspect will not be taken into account in this chapter. Even if active particles are obviously three-dimensional and their motion also happens in 3D, here we will mainly consider motion in two dimensions, as in many real situations the motion of active particles is in a quasi-2D environment; for example, the case of motile bacteria moving above the lower horizontal surface of a sample slide. However, we will provide details on how to handle the motion of active particles in 3D as well. It is worth noting that in some cases 2D-confinement on active particles can give rise Here are represented biological and artificial particles of micro-and nanoscopic size. Their speed ranges from a few µm s −1 to mm s −1 . Adapted from Bechinger et al. [2] to unexpected features that are not present in the case the active particles are moving in the bulk of a solution, i.e., full 3D motion with no confinement [3].

Passive Brownian Motion
Typical active particles are motile bacteria or artificial self-propelling microparticles, which perform their motion in a liquid environment. Since they are immersed in a fluid, active particles are subject to viscous force, always opposite to their velocity, and to thermal noise that is generated by the molecules constituting the fluid, which, because of the microscopic size of the particles, results in a non-negligible effect. Therefore, before entering into the details of how various models can describe the motion of active particles, we will discuss how the passive Brownian motion of a spherical particle can be described using Langevin equations and simulated using finite-difference equations [4].
Numerical simulations of Brownian dynamics date back to the 1970s and 1980s. Seminal works in this field are [5][6][7]. A comprehensive review is [8] and fundamental reference books are [9][10][11]. A reference for the numerical methods is [12], where the method with finite differences, better known as Euler-Maruyama scheme, is explained among other schemes. Here, we report the essential basics of passive Brownian dynamics and we invite the reader to refer to the milestone works cited above for a deeper insight and a comprehensive review.
Let us assume we have a spherical microscopic particle (for example, a transparent silica particle of 3 µm diameter) floating in a droplet of liquid solution (for example, a water solution) deposited on a microscopic glass slide. The mass of the particle is m ≈ 10 −14 kg. If we observe it with a microscope, we will see that the particle moves erratically, hovering above the flat glass surface of a microscope slide. If we track the particle by recording its position at times regularly separated by a fixed time interval t, we will find that its translational motion is purely diffusive, with translation diffusion constant for each of the two main directions D t given by where k B is the Boltzmann's constant, T is the absolute temperature, and γ t is the friction coefficient of the particle for translational displacements (in the bulk of a liquid solution, γ t = 6πηR, where η is the viscosity of the fluid and R the radius of the particle). This equation is the simplest expression of Einstein's fluctuationdissipation relation.
In the case of a homogeneous spherical particle with a perfectly smooth surface, it is not easy to experimentally detect the particle orientation. However, if we manage to measure the particle orientation, we will find that, in addition to an erratic translational motion, also the orientation of the particle changes randomly. If we can record the orientation, we will see that also the rotational motion is purely diffusive, this time with a different constant, the rotational diffusion constant D r given by where γ r = 8πηR 3 is the rotational friction coefficient of the particle. The cause of these translational and rotational erratic motions lies in the interactions of the suspended colloidal particle with the molecules constituting the fluids, which are affected by the temperature and, at equilibrium, present a velocity partition function distributed according to Maxwell's distribution [13]. Because of the collisions with the fluid molecules, the particle experiences a force and a torque that perturb its motion (thermal noise).
The translational dynamics of a particle in a fluid environment is described by the Langevin equation: where the term −γ t v is the viscous friction force of the fluid, and F th is the stochastic thermal force, which has zero average and variance 2 k B T γ t . Because of the tiny mass of a microscopic particle, often inertia can be neglected. In fact, the characteristic time needed to forget the inertial effects is the relaxation time τ rel = m/γ t , which increases with the mass of the particle m, and decreases with the friction coefficient γ t . For microscopic particles like the prototype silica particle of 2 µm diameter, the relaxation time τ rel is of the order of magnitude of 0.1 µs. Such relaxation time is several orders of magnitude below the time intervals typically sampled in experiments (for example, the time interval between two frames in an acquisition via a standard CMOS camera is of the order of some milliseconds). Therefore, Eq. (7.3) can be simplified to the overdamped Langevin equation: where the inertial term ma on the right side of Eq. (7.3) has been dropped. In all systems where t τ rel , Eq. (7.4) is sufficient to capture the relevant measurable physical features. In fact, it is possible to demonstrate that the solution of Eq. (7.3) converges to Eq. (7.4) in the limit m → 0 [14].
Often, Langevin equation (7.4) is rewritten as: dr = 2D t dW , (7.5) where dW is the derivative of a Wiener process, with zero average and variance 1 [4,15]. The simplest, though effective, way to a numerical solution of Eq. (7.5) is to use a finite-difference approach. In 2D, the variables in Eq. (7.5) can be explicitly written as: By writing the velocities as v x = x/ t and v y = y/ t, Eq. (7.6) takes the form [4,15] where W x and W y are realisations of independent stochastic processes with average 0 and standard deviation 1 [4,15]. By writing explicitly x = x n+1 − x n and y = y n+1 − y n , we are led to the finite-difference equation where, for a given t, we obtain the sequence {x n , y n } representing the trajectory of the passive Brownian particle in the plane.

Active Brownian Motion
One of the simplest models of active motion is active Brownian motion. Let us consider a spherical particle that self-propels with a constant speed v along a given internal orientation direction in 2D. Just like a passive Brownian particle, this particle is also affected by thermal noise, which affects both its translation and rotation. The configuration of the active Brownian particle is described by three variables: two spatial coordinates x and y for the position according to the lab reference frame, and one rotational coordinate θ for the orientation of the particle with respect to the lab reference frame [16]. The equations determining the dynamics are The finite-difference equations relative to Eqs. (7.9) are therefore the following: (7.10) Figure 7.2 depicts typical trajectories for an active Brownian particle with 3 µm diameter with characteristic self-propulsion speed between v = 0 µm s −1 (passive Brownian particle) and 12 µm s −1 . An important quantity for characterising the motion of microscopic systems is the mean square displacement (MSD). Considering a 2D motion, the MSD of a particle as a function of the time lapse t is where the average is often performed over time. 1 The MSD over a given time t represents, therefore, the average quadratic displacement from the position the particle had a time t before. From the MSD, we can gain a lot of insights about the dynamics of a system. The MSD of an active Brownian particle, i.e., for a trajectory governed by Eq. (7.9), is where t r = D −1 r is the characteristic time scale for the rotational diffusion, that, for a prototype particle with 3 µm diameter, is of about 20 s. If we write explicitly the expression above for the limits t t r and t t r , we have For a prototype particle with 3 µm diameter, the translational diffusion coefficient D t ≈ 0.3 µm 2 s −1 . If we consider an active particle with speed v = 5 µm s −1 , then v 2 = 25 µm 2 s −2 . If we consider the case t = 0.01t r = 0.2 s t r , then the diffusive contribution to the MSD (4 D t t ≈ 0.1 µm 2 ) is an order of magnitude smaller than the velocity contribution to the MSD (v 2 t 2 ≈ 1 µm 2 ). The two contributions become comparable, when t ≈ 0.001t r = 0.02 s. For smaller values of the time lapse t, the diffusive contribution prevails. However, if the reference active speed is faster than the v = 5 µm s −1 , as it happens in many cases (for example, compare typical speeds given in Fig. 7.1), and the size is bigger than 3 µm diameter (again, compare typical swimmers sizes given in Fig. 7.1), then the range for the time lapse where the velocity contribution prevails starts from few hundredths of second.
Therefore, the experimentally observed dependence on the time interval t is essentially ballistic (i.e., quadratic, ∝ v 2 t 2 ) for time scales smaller than the rotational diffusion time scale t r , and diffusive (i.e., linear, ∝ t) for time scales much longer than t r . In the latter case, the rotational diffusion plays a role in the randomisation of the propulsion direction over long times, and acts as an effective enhancement of the diffusion proportional to v 2 t r , so that the long-term effective diffusion coefficient is If there is no self-propulsion (i.e., v = 0), the dynamics is that of a passive Brownian particle and the MSD is linear at all time scales. Figure 7.3 represents the MSD corresponding to trajectories obtained from the numerical integration of Eq. (7.9), for different values of the self-propulsion velocity v.
It is an interesting exercise to calculate the MSD also for the cases of run-andtumble motion, chiral active Brownian motion, and active motion with Gaussian The MSD for the passive Brownian particle (v = 0 µm s −1 ) is also represented, and is linear over all the range of time delays noise reorientation mechanism, which will be described in the following subsections. Even though the microscopic locomotion mechanism differs at short time scales, the resulting MSDs present some universal features such as a ballistic behaviour at short time scales and an enhanced diffusive behaviour at long time scales.

Run-and-Tumble Motion
In the case of living organisms like motile bacteria, the observed motion can be described as a sequence of rectilinear forward steps and, occasionally, a sudden stop in the motion followed by a reorientation and by rectilinear motion along the new direction. Such kind of motion is known as run-and-tumble motion [1]. During the "runs", the bacterium moves forward because of the rotational motion of its flagella, which are tangled in a bunch. During the "tumbles", one of the flagella changes rotation direction breaking the bunch and the bacterium reorients itself. After this reorientation, the flagella form again a bunch and a new run starts. Run-andtumble motion is a typical strategy of chemotactic organisms, i.e., organisms that calibrate their motion according to the presence or absence of determined chemical substances, usually an attractant or a repellent, often following the concentration gradient. The mechanism of run-and-tumble motion has been thoroughly studied in E. coli [17,18] and various models have been developed to simulate it, from simple ones [19] to more complex ones [20].
In the absence of the chemical substance to which the chemotactic organism is responding, the reorientation events happen with a timing well-described by a Poisson process, characterised by the probability distribution: (7.15) where N is the number of events observed in the time interval t, λ is the average number of events expected in t, and n is a natural number. In this framework, the probability that tumbling happens is If we want to describe the system dynamics in a finite-difference equations approach, we have to add another variable to the set (x, y, θ), because we have to know if the bacterium is running or tumbling. Therefore, we add then the discrete variable that keeps track of the status of the bacterium (1: run, 0: tumble), and each time step t has a probability P tumble to be set to 0 and probability P run = 1−P tumble to be set to 1. The set of finite-difference equations is then: tumble,n n+1 = 0 or 1 with probability 1 − e −λ and e −λ each time step .
(7.17) An example of the appearance of a run-and-tumble trajectory is given in Fig. 7.4b.

Chiral Active Brownian Motion
It is not uncommon to observe that bacteria explore their environment by moving in circles. For example, E. coli bacteria have been shown to prefer to perform their quasi-circular motion in a clockwise fashion, when swimming close to a solid boundary, while they move counterclockwise when swimming near to an interface (for example, an air-liquid interface) [21][22][23][24][25]. The resulting motion is a chiral active Brownian motion. The term chiral, referring to the lack of symmetry of an object under mirror-reflection, is used to describe the tendency of swimmers to swim in circles either in a clockwise or a counterclockwise fashion. Therefore, the term chiral used in combination with the terms active Brownian motion indicates instead that the active particle tends to move in a well-defined way (clockwise or counterclockwise) with respect to a defined normal direction to the surface where the motion happens. Chiral motion is not only observed in living microorganisms, but also in artificial swimmers [26]. If we want to describe the dynamics of such a chiral active particle, we have to take into account that now the orientation varies with a defined angular velocity ω: According to the standard convention on the direction of angles, ω > 0 will be associated with a counterclockwise motion, while ω < 0 will characterise a clockwise motion. When translated to a finite-difference equation formalism, we have An example of the appearance of a chiral trajectory is given in Fig. 7.4c, d, for dextrogyre and levogyre chiralities, respectively.

Gaussian Noise Reorientation Model
Active motion is not only related to particles that are able to propel themselves in a strict sense. In fact, it has been observed that passive colloids in an active bath (for example, in a solution containing motile bacteria) present an effective dynamics that is quite different from standard passive Brownian dynamics. The presence in the solution of motile living microorganisms changes the motion of the suspended colloidal particles in such a way that they behave as effective active colloids for which the reorientation mechanism can present an enhanced diffusion constant [27][28][29][30]. The model commonly used for describing such a situation is This equation is practically Eq. (7.9), with the only difference that the noise term θ is not characterised by the rotational diffusion constant D r = k B T /(8πηR 3 ), but by a differentD r , usually larger than the one dictated by the size of the particle. Because the dependence of the noise term is still a Gaussian, the model with enhanced diffusion constant for the orientation is referred as Gaussian noise model. An example of the appearance of a Gaussian noise trajectory is given in Fig. 7.4e.

More Complex Models
Until now, we have learned to model the active motion of a single spherically symmetric particle in a 2D homogeneous environment. In this section, we will extend the model to include (1) the 3D case of a single spherically and nonspherically symmetric particle, (2) the presence of external fields, (3) the presence of multiple interacting particles, and (4) the presence of a multiplicative noise (that is, noise that depends on the state of the system and not only on external variables like the temperature). All these extensions are used widely in simulations for the description of the behaviour of real active systems. For example, one can have a colloidal solution of Janus particles that are concentrated enough to come close to each other, so that their (steric, electrostatic, hydrodynamic) interaction should be considered for a correct description of their collective behaviour [31][32][33] or a colloidal particle in a bacterial bath in the presence of an external optical potential inducing on the particles a driving force in the direction of the light intensity gradient [30]. Another common case are active dynamics of non-spherical particles like elongated rods [34][35][36][37] or chiral particles [38][39][40] in 3D or the presence of a boundary [41,42] that alters the value of the diffusion constant in its proximity via a diffusion gradient [43].

Non-Spherical Particles
If the particle is non-spherical, then the effect of the thermal noise is described by a diffusion matrix D of dimension 6 × 6, which takes into account all the translational and rotational modes of the non-spherical particle and possible correlations between them, including purely translational and purely rotational modes [44,45]. The diffusion matrix D is always symmetrical (D = D T ) and is represented as: where D tt is the diffusion term for the purely translational modes, D rr is the diffusion term for the purely rotational modes, and D tr and D rt are the off-diagonal terms describing roto-translational effect of the thermal agitation that might arise for particular particles shapes breaking mirror symmetry. With the same diffusion matrix formalism, we can describe the Brownian motion of a spherical particle, in which case D is a diagonal matrix, with D tt = D t I 3×3 and D rr = D r I 3×3 .
In case the passive motion of the spherical particle is in three dimensions, an equation analogous to Eq. (7.5) describes the dependence along all the three translation coordinates. However, because purely rotational and roto-translational terms may affect the orientation of the particle, it is a good practice to consider all the 6 degrees of freedom together. The analogue of Eq. (7.5) where the noise terms t , r are generated each step through a multi-varied Gaussian random number distribution satisfying the requirement of average equal to (0, 0, 0, 0, 0, 0) and variance matrix equal to 2 D t.
The increments r represent the displacement of the centre of mass of the particle with respect to the previous position.
The angular increments θ need to be handled with more care. The rotation of the particle is a free 3D rotation. Therefore, the increment θ represents a set of three increment, one for each rotation axis (that is, one for each of the unit vectors defining the set of axis of the particle reference frame) [45]. As the composition of rotations among different axes is not commutative, one should choose a time step t such that the various increments θ are small enough to ensure commutativity, within a certain error range. Each time step a rotation R = R x R y R z will be applied to the unit vectors determining the particle reference frame in order to find the rotated configuration. The rotation leaves the position of the centre of mass unaltered and, for the dynamics to make sense, should be such that the composition of the rotation matrices R x , R y , and R z is the same, within the order of t.
A cleaner approach to the issue takes into account the algebra of the generators of the rotation matrix group. A rotation of an angle φ around the x, y, z axis is written, in 3D, as a matrix acting on the component of a vector along the base unit vectors, and in the specific case is Each of this matrices can be written as an exponential of the generator matrices G x , G y , G z : in the following way: Instead of performing a small finite rotation about one axis at a time, where we might have issues because of non-commutativity, it is wiser to perform directly a rotation around the axis individuated by the vector ω of the angle θ = ( θ x ) 2 + ( θ y ) 2 + ( θ z ) 2 = t |ω| that is exactly the rotation acting on the particle. Such rotation matrix can be written as the exponential of the skew-symmetric matrix θ × : where Because of the properties of θ × , namely θ 3 × = −θ 2 θ × , can be written as: is the Rodrigues formula [46] for the rotation of an angle θ around a directionω that is exactly the rotation of the axes of the particle reference frame due to a rotational noise term (stochastic rotational torque) of ξ r .

External Fields
In many situations the particles feel the effects of external force or torque fields. In the case of colloidal particles suspended in a solution, these external fields can be due to the optical force generated by the presence of an optical potential [47], the presence of a hydrodynamic flux [48], the combined effect of weight and buoyancy that, for particles with an asymmetric mass distribution, can give rise to a torque leading to gravitaxis [49,50], or the presence on an external magnetic field for paramagnetic particles [51][52][53]. Moreover, in many current realisations of artificial microswimmers electric [54], magnetic [55,56], and acoustic fields [57,58], or a combination of them [59,60], play an important role to activate the selfpropelling mechanism, to control the swimming behaviour, or to confine the active particles [61]. A variety of models have been developed in order to understand the mechanism of the self-propulsion and its proper simulation [62][63][64][65]. When writing the equations of motion for the total force F tot and torque T tot acting on a particle, we need to include the contributions of the external fields F ext and T ext as: In the case of an active Brownian particle in 2D, the presence of an external force F ext and/or torque T ext,z is included in the equations as follows: In the general case of an active particle in a 3D environment, the effect of an external force F ext and torque T ext is included as follows: This formalism takes into account the effect of possible roto-translational effects, and reduces to the standard set of separable equations in the case translational and rotational motion are independent from each other. Besides the mathematical formalism, the presence of an external field might induce important features in the behaviour of a system of active Brownian particles. We will see an example of this in more details after discussing the modalities of simulating systems with more than one particle.

Interacting Particles
The study of the behaviour of a single active Brownian particle is an important starting point for understanding the behaviour of a system of multiple particles. When multiple particles are present, interactions among the particles may significantly change the dynamics. Such interactions, moreover, can affect the collective behaviour of the system, by determining the emergence of cooperative phenomena, like, for example, phase separation or the formation of dynamic clusters. Moreover, the presence of activity itself might change dramatically the behaviour of a system: the same interactions present in a system of active or passive particles may give rise to totally different outcomes, just because in one case we have the additional feature of self-propulsion. Usually interactions among active particles are divided into two main categories: non-aligning and aligning interactions. Non-aligning interactions might be attractive or repulsive, might depend on the relative position of the particles, but are independent on their direction of motion. Aligning interactions instead depend on the particles' direction of motion as well. Often such interactions tend to align the particles in their direction of motion, favouring phenomena such as swarming, like, for example, in the Vicsek model [66]. Theoretically, it has been predicted that active particles responding to chemical signalling or to hydrodynamic interactions may interact mutually with an effective aligning interaction [67][68][69][70].
We start our brief overview of interacting active particles with steric interactions, which prevent the particles from occupying the same volume, and we continue with two kind of aligning interaction, one characterising the Vicsek model, and another one giving a torque on particles at very short distance from one another.

Steric Interactions
Usually colloidal particles have a well-defined, rigid shape, and it is not possible for them to overlap. This steric interaction, that is present between passive and active colloids, together with the presence of active particles might give rise to interesting phenomena, like the formation of metastable clusters, even though the interaction itself is not attractive. For example, a set of passive colloids does not form spontaneously any cluster except that in the presence of a strong attractive interaction or of a driving force that pushes the passive particles all in the same space [30]. A set of active particles, instead, because of their activity can form metastable cluster in dilute suspensions even in the presence of repulsive mutual interactions [71][72][73][74]. In fact, depending on the propulsion velocity v and on the rotational diffusion time τ r , two active particles might collide together and stay locked because of the persistency of their active motion. Such clusters then can break apart over a time scale dictated by τ r because the reorientation process makes one of the particle to point away [75,76].
For what concerns the simulations of the cases where the particles are rigid and have a finite dimension, steric interactions can be implemented via the socalled hard-sphere correction, shown in Fig. 7.5, in order to avoid non-physical situations where the particles overlap and share in part the same volume. This is Fig. 7.5 Steric interactions. Schematic for hard-spheres interaction between two particles, where the particles are displaced away from each other whenever they get superposed done by checking the mutual distance between couples of particles after each time step: if an overlap is found (for example, when the distance between the particles' centres is smaller than the diameter of the particles), then both particles are moved away from each other along the direction connecting the centres in such a way that their distance becomes exactly one diameter, so they touch each other, but are not superposed. This procedure has been used in [77] and, although it does not implement an elastic collision among particles like the other approaches [9] where conservation of energy and momentum is enforced, it is computationally lighter. However, the allowed superposition among the superposing particles is usually a small fraction of the particle volume, so, the time step of the simulation should be accordingly small, not to lead to unphysical effects, like a state of excessive superposition or to the unphysical conditions of two nearby particles that, instead of colliding, miss each other because of an excessive value of the calculated displacement. In case numerous particles are in close contact to each other, like, for example, when they are part of a cluster, one should check iteratively each couple of neighbouring particles, so to ensure that, after each readjusting move, the condition of non-superposition remains valid. Such implementation has been shown to be equivalent to introduce a short-range repulsive potential that prevents superposition. However, computationally the hard-spheres correction is preferred, because for the correct dynamics with the short-range, repulsive potential one would have to use a much smaller time step, which is unpractical.
Vicsek Model The Vicsek model [66] is one of the simplest models to feature alignment and swarming in a system of active particles. In the Vicsek model, the particles move with constant speed and interact via the following aligning interaction: each particle can sense the orientation of the particles within a determined flocking radius R flocking , and at each step reorients itself according to the average value of the neighbouring particles' orientation ( Fig. 7.6). If, for each particle i in the system, we define S i as the set of neighbouring particles at the considered instant in time, then the equations describing the system are Varying the range of the parameters describing the system, within this model one can obtain a phase transition from undirected motion to unidirectional motion, as a function of the particles' density: beyond a given density, whatever the initial conditions on position and orientation, the particles will come to move in the same direction, thanks to the aligning interaction. If we include in Vicsek model the steric interaction, we can obtain crystallisation at high densities and low-noise intensity. Vicsek model: reorientation mechanism. The status at a given instant is represented. For establishing the orientation at the following instant of the particle represented in orange, one has to take the average of the orientations of the particles lying within the flocking radius (green area).
In the second panel the new orientation is shown in red. The previous orientation is shown in light blue, for a comparison

Short-Range Aligning Interactions
Here we present another mechanism of aligning interaction [78], which can describe the effective behaviour of particles interacting with some aligning hydrodynamic interactions, or of bacteria moving in a background of colloidal particles, or of people moving in a crowd. We have a set of finite-size active particles in 2D which, when closer than a given distance R align , interact by means of a torque. The behaviour featured by a particle is a tendency to reorient the direction of its active displacement towards the particles in the forward direction of motion, and away from the particles that are lying backwards. Ifk is the unit vector direction perpendicular to the plane of motion, the torque acting on a given particle will be alongk. The torque acting on a given particle n because of the presence of another particle i within the interaction distance is proportional to the cosine of the angle formed by the direction connecting the centre of the particle to the particle i (r ni ) with the direction of the velocity of the particle n (v n ), and inversely proportional to the square of the distance among the particles. In order to implement the behaviour described above, the torque should be along the vectorv n ×r ni , and proportional to (v n ·r ni )/r 2 ni . Such coefficient can be positive, negative, or zero, depending onv n ·r ni , and this is the ingredient that gives the wanted behaviour of reorienting towards particles set in the front, and turning away from the particles lying at the back (Fig. 7.7).
Expressing the total torque on a particle using the coefficient along the given k direction, we have that the total torque felt by the particle n can be modelled according to: Fig. 7.7 Short-range aligning interaction: reorientation mechanism. In (a) a finite set of relative positions and relative motions are shown. The particle represented in red (centre of the scheme) exerts a torque on the particle in black. Such torque obeys the mechanism in Eq. (7.36) and the reorientation angle it induces is represented, for the various case considered, by a red oriented arc. In (b-e) a few cases characteristic of the induced dynamics with two and three particles are represented. In particular, it is shown the situations leading to stable clusters of 2 particles (b, f) and 3 particles (g-i). In the cases represented in (c-e), the relative arrangement of the particles does not allow to form any cluster, but part away from each other (c, e) or proceed at constant speed and constant distance, with one of the particles following the other one. Reproduced from Nilsson and Volpe [78] (licenced under CC BY 3.0) T n = T 0 i∈Sv n ·r i r 2 iv ×r i ·k, (7.36) where S is the set of particles within the radius of interaction R align from particle n, and T 0 is a parameter setting the strength of the interaction. So, in general, the model is described by the following equations: ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ẋ n = v cos θ n + ξ x,ṅ y n = v sin θ n + ξ y,ṅ θ n = T n + ξ θ,n . (7.37) where T n is the torque obtained from Eq. (7.36) for each particle n. We note that in this specific model the rotational friction coefficient that usually is present in front of the torque is set to 1.
Such a model, though simple, displays a rich set of behaviours. In fact, depending on the rotational noise conditions and on the particle concentration, we can find a gaseous phase, where the particles are moving independently, a phase where metastable clusters composed by few particles are present, and a more pronounced clustering phase where clusters of more particles are formed. It has been shown [78] that the clustering transition occurs for the critical noise level of T 0 /4R 2 align . We expect also that the density should have an effect on the clustering transition; however, this aspect has not been investigated in the original paper. It has been also shown that in a mixed system of a few active and many passive particles, where the active particles interact with an alignment term as in Eq. (7.36) with the other active particles and with a term of the opposite sign where interacting with the passive particles, the interaction can lead to the formation of metastable channels [78].

Multiplicative Noise
In all the previous examples we have dealt with a uniform diffusion constant and therefore with noise conditions uniform in time, independent from the status of the system. However, in many real systems the noise may depend on the configuration of the system: for example, it is known that the presence of a rigid planar wall induces a gradient in the thermal diffusion coefficient, in such a way that the closer the particle to the planar boundary, the smaller the diffusion coefficient [79][80][81][82]. In particular, for the direction perpendicular to the planar surface, the diffusion coefficient D ⊥ becomes zero when the particle finds itself in contact with the wall [79,83]. An analogous effect happens also when two particles in the bulk of a solution come close to each other: the presence of one of the particle in the proximity of the other alters the diffusion coefficient of the second particle, and vice versa [84].
In general, a multiplicative noise is described as: ξ = (2 D(x))W, (7.38) where the diffusion constant depends on the value of the variable describing the status of the system. To fix our ideas, let us suppose it indicates the distance of a particle from the planar wall. In such cases, the full Langevin equation for the particle would be Eq. (7.3), with the only difference that the noise term has a dependence on the variable x. The corresponding equation is [15] where the additional term dD dx is the spurious drift and is necessary for the correct convergence of the solution to the original Langevin equation with multiplicative noise [14]. Finally, the corresponding finite-difference equation describing a particle in a gradient of diffusion is

Numerical Examples
Here we show two examples of collective behaviours emerging in a system of selfpropelling particles.

Living Crystals
We have performed a simulation of a system of N = 100 active particles considering steric interaction and a short-range phoretic attraction on the line of the model employed in [85], which was proposed to reproduce the experimentally observed living crystals emerging in a solution with light-activated colloids. In this model, the phoretic attractive force, due to the advective flow generated by the decomposition of hydrogen peroxide on the exposed hematite surface, induced an attractive velocity v P (r) for the active particles that scales according to the inverse of the square distance among the particles: Such dependence is the expected behaviour for a phoretic attraction to a reaction source, and fit well the observed experimental behaviour [85].
Varying the strength of this attractive interaction, which is the magnitude v P 0 of the attractive velocity at a given reference distance r 0 , and the self-propulsion velocity of the particle v, in the case particles are clustering, we can observe the formation of clusters of different average sizes. Such clusters are actively changing and rearranging themselves, and their stability depends also on the relative ratio of the interaction velocity at short distances (comparable with the diameter of the particles) with the self-propulsion velocity of the single active particle. For a given active velocity v, if the attractive reference velocity v P 0 for a given r 0 equal to the diameter of the particles is comparable to (or bigger than) the self-propulsion velocity v, then the clusters will tend to form. The bigger the reference value v P 0 of the attractive phoretic interaction, the more stable and bigger the cluster formed.
In Fig. 7.8 we provide two examples of simulations of active Brownian particles interacting via an attractive phoretic force that induces an attractive phoretic velocity given by Eq. (7.41). The radius of the particles is set to R = 1.0 µm, and their speed is set to v = 2 µm s −1 . In both the simulation presented, the initial configuration of positions and orientations is randomly generated. We consider a square arena with a side of L = 100 µm (concentration of φ V = 0.03) with periodic boundary conditions; the concentration of the particles can be adjusted by varying the size of the arena. The time step of the simulation is set to t = 0.1 s. With these parameters we observe the formation of metastable clusters as in Ref. [85]. By decreasing the strength of the interaction, we can see that the clusters become less stable. In Fig. 7.8 a few screenshots of the time evolution are presented.
In the dynamics represented in Fig. 7.8a-c, the reference phoretic velocity v P 0 is set to 2 µm s −1 at a distance between colloids of r 0 = 2 µm, and the radius of the Fig. 7.8 Living crystals. Simulation of 100 active particles of radius R = 1.0 µm with phoretic, short-range attractive interaction. The system forms metastable clusters. Depending on the particle concentration, on the noise conditions, and on the strength of the attractive interaction, the clusters may become more or less stable. In the simulation, the noise affecting the position and the orientation is set to the thermal noise felt by a spherical particle at equilibrium in a thermal bath. In (a-c), the phoretic interaction is set in such a way to give a phoretic attractive velocity equal to the speed of the active particle at a reference distance of 2.0 µm between the particles; the clusters formed are of the size of 10-20 particles, and are not too tightly bond and visibly evolve in the time scale of the simulation. In (d-f), the phoretic interaction was set to give a phoretic attractive velocity equal to the speed of the active particle at a reference distance of 2.5 µm between the particles; in this case the clusters formed are more compact, and, as the time pass, tend to form bigger clusters interaction is set to R int,phor = 10 µm [85]. We can observe that there is a tendency to form clusters, but they are not bigger than 20-30 particles each, and a consistent fraction of the particles are still freely moving. In the dynamics represented in Fig. 7.8d-f, instead, the reference phoretic velocity v P 0 is set to 2 µm s −1 at a reference distance between colloids of r 0 = 2.5 µm, i.e., the attractive interaction is stronger. Also in this case the radius of the interaction is set to R int,phor = 10 µm. As the attractive interaction is stronger, now the particles tend to form bigger clusters and less particles are freely moving independently.

Colloids with Short-Range Aligning Interaction
We consider a system of active particles interacting via the aligning interaction described in Eq. (7.36) and with a dynamics described by Eq. (7.37), following [78], where we explicitly take into account the proper value of the friction coefficient of a spherical particle in water for calculating the effect of the torque. We perform a simulation with N = 400 particles in a low-noise condition, with colloids of dimension R = 1.0 µm, setting the interaction radius R align = 2.5 µm, a speed of v = 0.1 µm s −1 , and torque T 0 = 0.5 rad µm 2 (Eq. (7.36)). With a simulation time step t = 0.05 s, we observe that clusters start forming and, because of the low-noise level, persist for a long time. In Fig. 7.9 a few screenshots of the time evolution are presented. In the chosen conditions of noise, range of the interaction, and strength of the interaction, we observe that mainly small clusters of two, three particles are formed. We observe also that some simple configurations like the ones shown in the circled regions of Fig. 7.9b-c (cluster of two particles facing each other with opposite orientation, and cluster of three particles organised in an equilateral triangle configuration with their orientation facing the centre of the triangle) are quite stable in time, and they do not move, and this happens because of As expected, small clusters of 2 and 3 particles are quite stable in time, given the low-noise level condition chosen, as can be seen comparing the configurations of the clusters evidenced in the circles the interaction mechanisms that tends to preserve the configuration and the relative low level of noise. Such small stable clusters break apart or transform only when a travelling particle or cluster comes in their proximity and, driven by the aligning interaction, targets them, as it is shown in the Supplementary Movies of Ref. [78]. Varying the conditions of noise, the strength of the interaction, the range of the interaction, one can reproduce various regimes, from a gaseous phase, with mainly small clusters and single isolated particles, to a cluster phase, with mainly clusters of various dimensions, the distribution of the size of the cluster being determined by the chosen condition of noise, strength of the interaction, and range of the interaction.

Conclusions
We have provided a concise, self-contained introduction to the simulation of active Brownian systems using Brownian dynamics and finite-difference equations. We have illustrated how to describe the dynamics of a micron-sized active particle starting from a model spherical particle, how to account for white thermal noise for translational and rotational motion, and active displacement, starting from a motion in two dimensions. We have briefly explained the way to proceed to simulate the full 3D motion, with the proper formalism also for non-spherically symmetric objects, and reviewed a few type of interactions that have been used in the description of the dynamics of active systems, and, depending on the choice of the parameters, present a variety of complex behaviours. We have tried to provide an overview that gives hints of the up-to-date research on active system, where simulation is often an important tool to test the understanding of the key mechanism determining the behaviour of the particular system under investigation. At the same time, we have kept the discussion as simple as possible and given a few numerical examples to provide an effective starting point for the students approaching numerical simulation of Brownian dynamics for the first time.