Lattice Boltzmann Methods and Active Fluids

TL;DR

This paper reviews the use of Lattice Boltzmann Methods in modeling active fluids, covering thermodynamics, implementation, and recent advances in phenomena like active turbulence and self-propelled droplets.

Contribution

It provides a comprehensive overview of how LBM can be applied to simulate complex active fluid systems and recent developments in the field.

Findings

LBM can accurately simulate active fluid hydrodynamics.

Recent studies have explored phenomena like active turbulence and self-propelled droplets.

Chapman-Enskog expansion links LBM to continuous active fluid equations.

Abstract

We review the state of the art of active fluids with particular attention to hydrodynamic continuous models and to the use of Lattice Boltzmann Methods (LBM) in this field. We present the thermodynamics of active fluids, in terms of liquid crystals modelling adapted to describe large-scale organization of active systems, as well as other effective phenomenological models. We discuss how LBM can be implemented to solve the hydrodynamics of active matter, starting from the case of a simple fluid, for which we explicitly recover the continuous equations by means of Chapman-Enskog expansion. Going beyond this simple case, we summarize how LBM can be used to treat complex and active fluids. We then review recent developments concerning some relevant topics in active matter that have been studied by means of LBM: spontaneous flow, self-propelled droplets, active emulsions, rheology, active turbulence, and active colloids.

Tables

Table 1. Table 1: The table summarizes bulk and elastic contributions to free energy for polar and nematic, both uniaxial and biaxial (see Appendix A), systems. Splay, twist and bending contributions have been written explicitly in terms of different elastic constants κ i subscript 𝜅 𝑖 \kappa_{i} ( i = 1 , 2 , 3 𝑖 1 2 3 i=1,2,3 ) for both polar and uniaxial nematic gels, while in the most general case of a biaxial nematic we did not distinguish between different contributions. The last line in the Table shows how the elastic contribution looks like assuming that the medium is elastically isotropic, i.e. , κ 1 = κ 2 = κ 3 = κ subscript 𝜅 1 subscript 𝜅 2 subscript 𝜅 3 𝜅 \kappa_{1}=\kappa_{2}=\kappa_{3}=\kappa .

Free energy contributions		Polar Gel	Nematic Gel
			Uniaxial	Biaxial
Bulk		$a{\overrightarrow{P}}^2+b{\overrightarrow{P}}^4$	$r{S}^2-w{S}^3+u{S}^4$	$\tilde{r}{Q}_{ij}{Q}_{ji}-\tilde{w}{Q}_{ij}{Q}_{jk}{Q}_{ki}+\tilde{u}{\left(Q_{ij}{Q}_{ji}\right)}^2$
Elastic	Splay	$\frac{{\kappa }_1}{2}{\left(\nabla \cdot \overrightarrow{P}\right)}^2$	$\frac{{\kappa }_1}{2}{\left(\nabla \cdot \overrightarrow{n}\right)}^2$
	Twist	$\frac{{\kappa }_2}{2}{\left(\overrightarrow{P}\cdot \nabla \times \overrightarrow{P}\right)}^2$	$\frac{{\kappa }_2}{2}{\left(\overrightarrow{n}\cdot \nabla \times \overrightarrow{n}\right)}^2$	$\frac{L_1}{2}{\left({\partial }_k{Q}_{ij}\right)}^2+\frac{L_2}{2}{\left({\partial }_j{Q}_{ij}\right)}^2+\frac{L_3}{2}{Q}_{ij}\left({\partial }_i{Q}_{kl}\right)\left({\partial }_j{Q}_{kl}\right)$
	Bend	$\frac{{\kappa }_3}{2}{\left(\overrightarrow{P}\times \nabla \times \overrightarrow{P}\right)}^2$	$\frac{{\kappa }_3}{2}{\left(\overrightarrow{n}\times \nabla \times \overrightarrow{n}\right)}^2$
Single constant approximation		$\kappa {\left(\nabla \overrightarrow{P}\right)}^2$	$\kappa {\left(\nabla \overrightarrow{n}\right)}^2$	$L_1{\left({\partial }_k{Q}_{ij}\right)}^2$

Table 2. Table 2: Explicit expressions of the elastic (first row) and the interface (second row) stress, and of the term S → → 𝑆 \vec{S} in the Beris-Edwards equation ( 18 ) (fourth row) for polar and nematic gels. The molecular field 𝚵 𝚵 \mathbf{\Xi} is a vector, with components h α subscript ℎ 𝛼 h_{\alpha} , for polar gels and a tensor H α ​ β subscript 𝐻 𝛼 𝛽 H_{\alpha\beta} , for nematic gels, as shown in the third row. κ 𝜅 \kappa is the elastic constant of the liquid crystal; the flow-alignment parameters ξ 𝜉 \xi and ξ ′ superscript 𝜉 ′ \xi^{\prime} are respectively related to the polarization field P → → 𝑃 \vec{P} and to the nematic tensor Q ¯ → ¯ ¯ ¯ → 𝑄 \underline{\underline{\vec{Q}}} and depend on the geometry of the microscopic constituents (for instance ξ > 0 𝜉 0 \xi>0 , ξ < 0 𝜉 0 \xi<0 and ξ = 0 𝜉 0 \xi=0 for rod-like, disk-like and spherical particles, respectively). In addition, these parameters establish whether the fluid is flow aligning ( | ξ | > 1 𝜉 1 |\xi|>1 ) or flow tumbling ( | ξ | < 1 𝜉 1 |\xi|<1 ) under shear. D ¯ ¯ = ( W ¯ ¯ + W ¯ ¯ T ) / 2 ¯ ¯ 𝐷 ¯ ¯ 𝑊 superscript ¯ ¯ 𝑊 𝑇 2 \underline{\underline{D}}=(\underline{\underline{W}}+\underline{\underline{W}}^{T})/2 and Ω ¯ ¯ = ( W ¯ ¯ − W ¯ ¯ T ) / 2 ¯ ¯ Ω ¯ ¯ 𝑊 superscript ¯ ¯ 𝑊 𝑇 2 \underline{\underline{\Omega}}=(\underline{\underline{W}}-\underline{\underline{W}}^{T})/2 represent the symmetric and the antisymmetric parts of the velocity gradient tensor W α ​ β = ∂ β v α subscript 𝑊 𝛼 𝛽 subscript 𝛽 subscript 𝑣 𝛼 W_{\alpha\beta}=\partial_{\beta}v_{\alpha} .

	Polar Gel	Nematic Gel
${\sigma }_{\alpha \beta }^{elastic}$	$\frac{1}{2}\left(P_{\alpha }{h}_{\beta }-P_{\beta }{h}_{\alpha }\right)-\frac{\xi }{2}\left(P_{\alpha }{h}_{\beta }+P_{\beta }{h}_{\alpha }\right)-\kappa {\partial }_{\alpha }{P}_{\gamma }{\partial }_{\beta }{P}_{\gamma }$	$2{\xi }^{\prime }\left(Q_{\alpha \beta }-\frac{{\delta }_{\alpha \beta }}{3}\right)Q_{\gamma \nu }{H}_{\gamma \nu }-{\xi }^{\prime }{H}_{\alpha \gamma }\left(Q_{\gamma \beta }+\frac{{\delta }_{\gamma \beta }}{3}\right)$ $-{\xi }^{\prime }\left(Q_{\alpha \gamma }+\frac{{\delta }_{\alpha \gamma }}{3}\right)H_{\gamma \beta }-{\partial }_{\alpha }{Q}_{\gamma \nu }\frac{\delta ℱ}{\delta {\partial }_{\beta }{Q}_{\nu \gamma }}+Q_{\alpha \gamma }{H}_{\gamma \beta }-H_{\alpha \gamma }{Q}_{\gamma \beta }$
${\sigma }_{\alpha \beta }^{interface}$	$\left(f-\varphi \frac{\delta F}{\delta \varphi }\right){\delta }_{\alpha \beta }-\frac{\partial f}{\partial \left({\partial }_{\beta }\varphi \right)}{\partial }_{\alpha }\varphi$	$\left(f-\varphi \frac{\delta F}{\delta \varphi }\right){\delta }_{\alpha \beta }-\frac{\partial f}{\partial \left({\partial }_{\beta }\varphi \right)}{\partial }_{\alpha }\varphi$
$𝚵$	$h_{\alpha }=\frac{\delta F}{\delta P_{\alpha }}$	$H_{\alpha \beta }=\frac{\delta F}{\delta Q_{\alpha \beta }}-\left(\frac{\delta F}{\delta Q_{\gamma \gamma }}\right){\delta }_{\alpha \beta }$
$\overrightarrow{S}$	$-{\Omega }_{\alpha \beta }{P}_{\beta }+\xi D_{\alpha \beta }{P}_{\beta }$	$\left[{\xi }^{\prime }{D}_{\alpha \gamma }+{\Omega }_{\alpha \gamma }\right]\left(Q_{\gamma \beta }+\frac{{\delta }_{\gamma \beta }}{3}\right)+\left(Q_{\alpha \gamma }+\frac{{\delta }_{\alpha \gamma }}{3}\right)\left[{\xi }^{\prime }{D}_{\gamma \beta }-{\Omega }_{\gamma \beta }\right]$ $-2{\xi }^{\prime }\left(Q_{\alpha \beta }+\frac{{\delta }_{\alpha \beta }}{3}\right)\left(Q_{\gamma \delta }{\partial }_{\gamma }{v}_{\delta }\right)$

Table 3. Table 3: Lattice speeds with their weights ω i subscript 𝜔 𝑖 \omega_{i} for spatial dimensions d = 2 𝑑 2 d=2 and d = 3 𝑑 3 d=3 and number of neighboring nodes Q 𝑄 Q .

Lattice	${\overrightarrow{\xi }}_i$	${\omega }_i$
$d2Q7$	$(0,0)$	$1/2$
	$c(\cos \left(i\pi /3\right),\sin \left(i\pi /3\right))$	$1/12$
$d2Q9$	$(0,0)$	$4/9$
	$(\pm c,0)(0,\pm c)$	$1/9$
	$(\pm c,\pm c)$	$1/36$
$d3Q15$	$(0,0,0)$	$2/9$
	$(\pm c,0,0)(0,\pm c,0)(0,0,\pm c)$	$1/9$
	$(\pm c,\pm c,\pm c)$	$1/72$
$d3Q19$	$(0,0,0)$	$1/3$
	$(\pm c,0,0)(0,\pm c,0)(0,0,\pm c)$	$1/18$
	$(\pm c,\pm c,0)(\pm c,0,\pm c)(0,\pm c,\pm c)$	$1/36$
$d3Q27$	$(0,0,0)$	$8/27$
	$(\pm c,0,0)(0,\pm c,0)(0,0,\pm c)$	$2/27$
	$(\pm c,\pm c,0)(\pm c,0,\pm c)(0,\pm c,\pm c)$	$1/54$
	$(\pm c,\pm c,\pm c)$	$1/216$

Equations

$$\vec{P}(\vec{r},t)=⟨\vec{\nu }(\vec{r},t)⟩=\int \mathrm{d\Omega }{f}_P(\vec{\nu },\vec{r},t)\vec{\nu },$$

$$\vec{P}(\vec{r},t)=P(\vec{r},t)\vec{n}(\vec{r},t),$$

$$Q_{\alpha \beta }(\vec{r},t)=⟨{\nu }_{\alpha }{\nu }_{\beta }-\frac{1}{d}{\delta }_{\alpha \beta }⟩=\int \mathrm{d\Omega }{f}_Q(\vec{\nu },\vec{r},t)({\nu }_{\alpha }{\nu }_{\beta }-\frac{1}{d}{\delta }_{\alpha \beta }).$$

$$Q_{\alpha \beta }(\vec{r},t)=S(\vec{r},t)\left[n_{\alpha }(\vec{r},t)n_{\beta }(\vec{r},t)-\frac{1}{d}{\delta }_{\alpha \beta }\right].$$

$$S(\vec{r},t)=\frac{1}{2}⟨3{\cos }^2\theta -1⟩,$$

$$\begin{array}[]{l}\displaystyle L_{1}=\dfrac{\kappa_{3}+2\kappa_{2}-\kappa_{1}}{9S^{2}},\\ \\ \displaystyle L_{2}=\dfrac{4(\kappa_{1}-\kappa_{2})}{9S^{2}},\\ \\ \displaystyle L_{3}=\dfrac{2(\kappa_{3}-\kappa_{1})}{9S^{3}},\end{array}$$

$${\sigma }_{\alpha \beta }^{\text{active}}=-\zeta \varphi Q_{\alpha \beta },$$

$${\sigma }_{\alpha \beta }^{\text{active}}=-\zeta \varphi \left(P_{\alpha }{P}_{\beta }-\frac{1}{d}|\vec{P}{|}^2{\delta }_{\alpha \beta }\right).$$

$$\underline{\underline{\sigma }}=({\Gamma }_2{\nabla }^2+{\Gamma }_4{\nabla }^4)\left[\nabla \vec{v}+(\nabla \vec{v}{)}^T\right].$$

$$F$$

$$\mu =\frac{\delta F}{\delta \varphi }=a\varphi +b{\varphi }^3-k{\nabla }^2\varphi ,$$

$${\sigma }_{\alpha \beta }^a=-\hat{\zeta }\left({\partial }_{\alpha }\varphi {\partial }_{\beta }\varphi -\frac{1}{d}{\left(\nabla \varphi \right)}^2{\delta }_{\alpha \beta }\right),$$

$${\partial }_t\rho +\nabla \cdot (\rho \vec{v})$$

$$\rho \left({\partial }_t+\vec{v}\cdot \nabla \right)\vec{v}$$

$$\nabla \cdot \vec{v}=0+\mathcal{O}(M{a}^2),$$

$$\underline{\underline{\sigma }}={\underline{\underline{\sigma }}}^{\text{passive}}+{\underline{\underline{\sigma }}}^{\text{active}}.$$

$${\underline{\underline{\sigma }}}^{\text{passive}}={\underline{\underline{\sigma }}}^{\text{viscous}}+{\underline{\underline{\sigma }}}^{\text{elastic}}+{\underline{\underline{\sigma }}}^{\text{interface}}\text{mbox}.$$

$${\sigma }_{\alpha \beta }^{viscous}=\eta ({\partial }_{\alpha }{v}_{\beta }+{\partial }_{\beta }{v}_{\alpha })+\left(\tilde{\zeta }-\frac{2\eta }{d}\right){\partial }_{\gamma }{v}_{\gamma }{\delta }_{\alpha \beta },$$

$$\left({\partial }_t+\vec{v}\cdot \nabla \right)\mathbf{\Psi }-\vec{S}=-\Gamma \mathbf{\Xi },$$

$$D_t\overrightarrow{\stackrel{~}{r}}={\partial }_t\overrightarrow{\stackrel{~}{r}}+(\vec{v}\cdot \nabla )\overrightarrow{\stackrel{~}{r}}+\underline{\underline{D}}\cdot \overrightarrow{\stackrel{~}{r}}+\underline{\underline{\Omega }}\cdot \overrightarrow{\stackrel{~}{r}},$$

$${\partial }_t\varphi +\nabla \cdot \left(\varphi \vec{v}\right)=\nabla \cdot \left(M\nabla \frac{\delta F}{\delta \varphi }\right),$$

$$f_i(\vec{r}+{\vec{\xi }}_i\Delta t,t+\Delta t)-f_i(\vec{r},t)=\mathcal{C}(\{f_i\},t),$$

$${\partial }_t{f}_{ki}+({\vec{\xi }}_{ki}\cdot \nabla )f_{ki}=\mathcal{C}(\{f_{ki}\},t).$$

$$\mathcal{C}(\{f_i\},t)=-\frac{1}{\tau }(f_i-f_i^{eq}),$$

$$\rho (\vec{r},t)$$

$$\rho (\vec{r},t)\vec{v}(\vec{r},t)$$

$$\rho (\vec{r},t)$$

$$\rho (\vec{r},t)\vec{v}(\vec{r},t)$$

$$f_i^{eq}=A_s+B_s{v}_{\alpha }{\xi }_{i\alpha }+C_s{v}^2+D_s{v}_{\alpha }{v}_{\beta }{\xi }_{i\alpha }{\xi }_{i\beta },$$

$$\sum _i{f}_i^{eq}{\xi }_{i\alpha }{\xi }_{i\beta }=\frac{c^2}{3}\rho {\delta }_{\alpha \beta }+\rho v_{\alpha }{v}_{\beta }.$$

$$A_0=\rho -20A_2$$

$$B_0=0$$

Full text

1 Introduction

The goal of this article is to describe the use of Lattice Boltzmann Methods in the study of large scale properties of active fluids toner2005 ; ramaswamy2010 ; marchetti2013 ; elgeti2015 ; gonnella2015 ; bechinger2016 ; doostmohammadi2018 , also showing recent progress in few relevant topics. Active fluids are living matter or biologically inspired systems with the common characteristic of being composed by self-propelled (or active) units that burn stored or ambient energy and turn it into work giving rise, eventually, to systematic movement. An example in nature is given by the cell cytoskeleton or, in laboratory, by synthetic suspensions of cell extracts with molecular motors (e.g. myosin or kinesin) surrey2001 ; bendix2008 . Molecular motors exert forces on cytoskeletal filaments (actin filaments and microtubules) howard2003 and trigger their motion in the surrounding fluid. These forces, exchanged through transient and motile contact points between filaments and motor proteins, result from the conversion of chemical energy, typically coming from ATP hydrolysis, into mechanical work.

Active systems show many interesting physical properties, of general character, related to their collective behavior, remarkable especially when compared with their analogue in passive or equilibrium systems. Pattern formation is an example. A disordered array of microtubules may arrange into spiral or aster configurations when the concentration of motor proteins like kinesin is sufficiently high surrey2001 . Suspensions of bacteria, despite their low Reynolds numbers, can exhibit turbulent flow patterns dombrowski2004 ; dunkel2013 , characterized by traveling jets of high collective velocities and surrounding vortices. Active fluids can be classified according to their swimming mechanism as extensile or contractile, if they respectively push or pull the surrounding fluid. This difference marks all the phenomenology of active fluids and, in particular, has important effects on the rheological properties. Activity is either capable to develop shear-thickening properties in contractile systems hatwalne2004 ; marenduzzo2007 ; foffano2012 ; foffanoPRL ; sokolov2009 , or to induce a superfluidic regime under suitable conditions in extensile suspensions cates2008 ; loisy2018 ; Guo201722505 . Simulations of extensile active emulsions under constant shear have shown the occurrence of velocity profiles (for the component of velocity in the direction of the applied flow) with inverted gradient (negative viscosity) and also jumps in the sign of apparent viscosity cates2008 ; loisy2018 ; Guo201722505 .

Other striking properties have emerged in the study of fluctuation statistics simha2002 ; chate2006 ; narayan2007 ; deseigne2010 ; kumar2011 ; cagnetta2017 and of order-disorder phase transitions  bechinger2016 ; cates ; paoluzzi2016 . Fluctuations and phase transitions have been mainly analyzed in the context of agent-based models. The flocking transition cavagna2014 , for instance, was the first one to be studied in a model of point-like particles moving at fixed speed and with aligning interaction vicsek1995 . Activity alone actually favors aggregation and can induce a phase transition, often called Motility Induced Phase Separation (MIPS) cates2015 . This has been numerically studied by using simple models of active colloids with excluded volume interactions and various shapes peruani2006 ; baskaran2008 ; yang2010 ; fily2012 ; lowen2013 ; redner2013 ; cugliandolo2017 ; digregorio2018 . The particle description has been also largely used in other contexts, to simulate, for example, the self-organization of cytoskeleton filaments described as semiflexible filaments nedelec1997 .

By a different approach, large scale behavior and macroscopic material properties of active fluids have been largely studied using coarse-grained descriptions based on general symmetry arguments and conservation laws. The first continuum description in terms of density and polarization field, with interactions favoring alignment with polar order, was proposed in toner1995 . In this model, as in others where nematic interactions were considered gruler1999 ; ramaswamy2003 ; baskaran2010 , the medium in which particles are supposed to move does not contribute with its own dynamics to the evolution of the system. Hence the environment of the active system can be considered as a momentum-absorbing substrate so that momentum is not conserved. On the other hand, there are systems in which the dynamics of the solvent can be relevant in a certain interval of length scales schaller2010 and must be incorporated in the description. The action of the active components on the solvent is taken into account by introducing an active stress into a generalized form of the Navier-Stokes equation. Suitable advection terms depending on the self-propulsion velocity of active units also appear in the dynamical equations for the order parameters describing the orientation of the active material (nematic or polar) or its concentration. A useful form for the active stress was first proposed in pedley1990 and later developed in the context of a coarse-grained model in simha2002 and, for active filaments or orientable particles, in kruse2004 ; julicher2007 . The total stress also includes elastic contributions, depending on the polar or nematic character of the system, stemming from an appropriate free-energy expression, as in the passive or equilibrium counterpart of the systems in exam usually called active gels. The resulting dynamical description consists of non-linear coupled partial differential equations that require numerical methods to be solved. A suitable approach, largely used to study multicomponents and complex fluids whose dynamics obey such equations, is the Lattice Boltzmann Method (in the following we will refer to it as LBM or LB) higuera1989 ; succi1991 ; succi2001 , a computational fluid dynamics scheme for solving the Navier-Stokes equation, eventually coupled to advection-relaxation equations yeomans2009 ; swift1996 ; gonnella1997 ; lamura1999 ; tiribocchi2009 ; denniston2001 . Among its features, this method is found to correctly capture the coupling between hydrodynamics and orientational order of liquid crystals (often known as backflow leslie1968 ; martin1972 ; beris1994 ), a crucial requirement to simulate the dynamics of active gels orlandini2008 ; marenduzzo2007 .

In this article we will review the way LBM can be used to describe collective properties of active fluids, describing also recent developments concerning issues where hydrodynamics plays a relevant role. We will initially review the thermodynamics of active fluids whose internal constituents are orientable objects, such as active liquid crystals. After shortly introducing the order parameters and the free-energy usually adopted to describe their properties, we will show how the active behavior enters the model and how hydrodynamic equations can be written to correctly capture the physics. This will be done in Section 2.

Afterwards we will discuss different LB strategies used to study simple and structured fluids, convenient for active fluids generalization. For a simple fluid, LBM solves a minimal Boltzmann kinetic equation governing the evolution of a single set of variables (the distribution functions), in terms of which hydrodynamic quantities can be written mcNamara1988 ; higuera1989 . A detailed description of the LB methods for a single fluid can be found in succi2001 ; Chen1992 ; succi2018 . For structured fluids, a full LBM approach can be followed by introducing a further set of distribution functions for the order parameter that follow the dynamics of appropriate lattice Boltzmann equations to be added to those describing the dynamics of the density and velocity of the fluid swift1996 . Then, interactions can be implemented by specific collision rules introduced on a phenomenologically ground or by making reference to a specific free-energy model that sets the thermodynamics of the system swift1995 ; swift1996 ; li2012 ; li2012_2 . The first approach, in numerous variants, has been largely used in the context of binary mixtures, due to its practical convenience, with the collision step designed in order to favor separation of the A and B components of the mixture shan1993 . When the fluid structure becomes more complex, the second approach becomes almost mandatory. The characteristics of a specific system will enter the lattice dynamic equations through a chemical potential and a pressure tensor that can be obtained by a given free-energy functional. Liquid crystals denniston2001 , but also ternary mixtures with surfactant lamura1999 or other kinds of complex fluids Cates2099 ; tiribocchi2011 , have been largely studied in this way. Finite difference methods, with possible numerical advantages, can be also applied to simulate the order parameter dynamical equations tiribocchi2009 and have been implemented in hybrid approaches coupled to LBM used as a solver for Navier-Stokes equations. These different options will be reviewed in Section 3, in relation with the modeling of active fluids proposed in Section 2, and with details on possible algorithms and numerical implementations.

The following Sections will be dedicated to discuss some relevant topics in active fluids in which LBM has played an essential role. In Section 4 the main numerical results concerning the hydrodynamic instabilities generated by spontaneous flows voituriez2005 ; marenduzzo2010 will be reviewed. Understanding how this occurs is fundamental, for instance, to assess the dynamics of topological defects as well as the physics of self-propelled droplets, objects which can capture some relevant features of motile cells tjhung2012 ; tjhung2015 . Section 5 will be devoted to review relevant results on the modeling of self-propelled droplets and of systems with many droplets such as active emulsions. The latter is a new subject of research with new fascinating perspectives. Active emulsions can be potentially realized by dispersing sticky bacteria linek2012 or self attractive cytoskeleton gels sanchez2012 ; guillamat2016 in water, or encapsulating an active nematic gel within a water-in-oil emulsion sanchez2012 ; guillamat2016 . Another stimulating field of research concerns the study of the rheological response of an active fluid to externally imposed flows. In Section 6 we will review the most recent and pioneering achievements in this field, in which, for example, an active gel has been predicted to have either a shear-tickening or superfluid-like behavior depending on the nature, extensile or contractile, of the flow PhysRevE.83.041910 . As a further topic we will illustrate the results obtained via LBM simuations to investigate “active” turbulence dombrowski2004 ; wensink2012 ; creppy2015 ; putzig2016 , a turbulent-like behavior observed in active fluids at low Reynolds numbers (Section 7).

The versatility of LBM as a solver for hydrodynamics has been also used to study the flow generated by different kinds of swimmers, treated as discrete particles coupled to the surrounding fluid by proper boundary conditions pooley2008 . In this case the solvent is described as a simple fluid whose dynamics can be solved by LBM. Although the study of the collective properties of these systems can be difficult within this approach due to the complicated structure of the flows induced by the swimmers, in a few cases this shortcoming has been overcome by using a mixed particle-continuum description ramachandran2006 ; Llopis2006 ; alarcon2013 ; earl2007 ; smith2007 ; degraaf2016 ; degraaf2016_2 . Section 8 will be dedicated to describe how LBM has been extended to include active particles.

2 Active fluid models

In this section we will focus on fluids whose internal units have an orientable character, a feature that crucially affects their reciprocal interactions, especially when a high density sample of active units is considered. In such cases the emerging orientational order on macroscopic scales can be captured by proper order parameters, such as the polarization vector $\vec{P}(\vec{r},t)$ and the tensor ${\underline{\underline{\vec{Q}}}}(\vec{r},t)$, often used to describe ordering in liquid crystals. These quantities will be introduced in 2.1.

The thermodynamics of these systems is usually described via a Landau-like free energy functional, depending upon powers of the order parameter and its gradients, respecting the symmetries of the disordered phase. The different free-energy terms describing bulk and elastic properties of the active fluid will be discussed in 2.2, while 2.3 will be dedicated to describe how activity is introduced in continuum models. In 2.4 we will briefly discuss the thermodynamics of a fluid mixture with an active component, with and without alignment interaction. The latter case has been recently considered for the study of the motility induced phase separation in active fluids tiribocchi2015 .

Finally the hydrodynamic equations describing both the evolution of the order parameter and of the velocity field will be shown in 2.5.

2.1 Order parameters

Active fluids whose internal constituents have an anisotropic shape (such as an elongated structure) encompass diverse systems ranging from bacterial colonies and algae suspensions marchetti2013 to the cytoskeleton of eukaryotic cells bray2000 . Depending upon the symmetries of such microscopic agents and upon their reciprocal interactions, these active fluids generally fall into two wide categories. The first one is the active polar fluid composed of elongated self-propelled particles, characterized by a head and a tail, whose interactions have polar symmetry. Such systems may order either in polar states, when all the particles are on average aligned along the same direction, as in the case of bacteria self-propelled along the direction of their head dombrowski2004 . Nevertheless systems of intrinsic polar particles, such as actin filaments cross-linked with myosin bray2000 ; tjhung2011 ; tjhung2012 ; nature2015 or microtubule bundles coupled with kinesin motors surrey2001 ; Decamp2015 ; sanchez2012 , may still arrange in a nematic fashion, restoring head-tail symmetry, when interactions favour alignment regardless of the polarity of the individual particles. Fig. 1 shows, for example, the aggregated phase of a system of self-propelled Brownian polar dumbbells suma2014 ; suma2014_2 ; gonnella2014 ; cugliandolo2015 which, depending on the strength of the self-propulsion force, may arrange in a polar state (right) or in an isotropic state (left), a behavior also found in bacterial colonies dellarciprete2018 . The second class includes head-tail symmetric, or apolar, particles that may move back and forth with no net motion, and order in nematic states. Examples of realizations in nature include melanocytes Kemkemer2000 , i.e. melanin producing cells in human body, and fibroblasts Duclos2014 , cells playing a central role in wound healing, both spindle-shaped with no head-tail distinction.

The continuum fields describing polar and nematic order are the vector field $P_{\alpha}(\vec{r},t)$ and the tensor field $Q_{\alpha\beta}(\vec{r},t)$ respectively (Greek subscripts denote the Cartesian components). They emerge either from a coarse grained description of a microscopical model degennes1993 or from a theory based on general symmetry arguments martin1972 ; beris1994 . Following, for instance, the former approach, for a system of rod-like particles the polarization field can be defined as

[TABLE]

where $f_{P}(\vec{\nu},\vec{r},t)$ is the probability density, encoding all the information coming from the microscopical model, of finding a particle at position $\vec{r}$ and at time $t$ oriented along the direction $\vec{\nu}$, and the integration is carried out over the solid angle $\mathrm{\Omega}$. The polarization can be also written as

[TABLE]

where $\vec{n}(\vec{r},t)$ is a unit vector defining the local mean orientation of particles in the neighborhood of $\vec{r}$, and $P(\vec{r},t)$ is a measure of the local degree of alignment, ranging from [math] (in an isotropic state) to $1$ (in a perfectly polarized state).

Differently, the nematic phase cannot be described by a vector field, as both orientations $\vec{\nu}$ and $-\vec{\nu}$ equally contribute to the same ordered state, due to the head-tail symmetry of the constituents. For a system of rod-like particles, the order is described by a nematic tensor which, in the uniaxial approximation (i.e. when a liquid crystal is rotationally symmetric around a single preferred axis), can be defined as

[TABLE]

Again $f_{Q}(\vec{\nu},\vec{r},t)$ is the probability density to find a nematic particle oriented along $\vec{\nu}$ at position $\vec{r}$ and time $t$, while $d$ is the dimensionality of the system. As for the polarization field, the nematic tensor can be also written in terms of the versor $\vec{n}$ (usually called director field) defining the local mean orientation of the particles

[TABLE]

Note that, by defining the nematic tensor in this way, one can separate local anisotropic features out of isotropic ones. Indeed, the only scalar quantity that can be derived from a tensorial object, i.e. its trace, is identically null. In Eq. (4) $S(\vec{r},t)$ plays the same role of $P(\vec{r},t)$ in defining the degree of alignment of the molecules in the nematic phase. In fact, by multiplying Eq. (3) and (4) by $n_{\alpha}n_{\beta}$, summing over spatial components and comparing them, one gets (in three dimensions)

[TABLE]

where $\cos\theta=\vec{n}\cdot\vec{\nu}$ is a measure of the local alignment of particles. The scalar order parameter $S$ achieves its maximum in the perfectly aligned state, where $\langle\cos^{2}\theta\rangle=1$, while it falls to zero in the isotropic phase where the probability density $f_{Q}$ is uniform over the solid angle and $\langle\cos^{2}\theta\rangle=1/3$. Assuming $\vec{n}$ to be parallel to a Cartesian axis, one can soon verify from Eq. (4) that $Q_{\alpha\beta}$ has two degenerate eigenvalues $\lambda_{2}=\lambda_{3}=-S/3$ (whose associated eigenvectors lie in the plane normal to the particle axes) and a third non-degenerate one $\lambda_{1}=2S/3$, greater in module than $\lambda_{2}$ and $\lambda_{3}$ and related to the director itself. Such formalism can be also extended to treat the case of biaxial nematics, i.e. liquid crystals with three distinct optical axis. Unlike an uniaxial liquid crystal which has an axis of rotational symmetry (such as the director ${\vec{n}}$), a biaxial liquid crystal has no axis of complete rotational symmetry. As such theory is out of the scope of this review, we briefly mention it in Appendix A, focusing, in particular, on how biaxiality is included in the tensor order parameter and on the role it plays in the localization of topological defects.

2.1.1 Topological defects

Topological defects (disclinations in nematic/polar and cholesteric liquid crystals) are regions where the order parameter cannot be defined degennes1993 ; chaikin1995 . A crucial difference between allowed polar and nematic systems really lies on the nature of the topological defects. As they play a relevant role in the dynamics of the velocity field in active fluids, we provide here a brief introduction about the theory of topological defects and address the reader to more specialized books (such as chaikin1995 ) for further details.

A topological defect can be characterized by looking at the configuration of the order parameter far from its core. This can be done by computing the winding number (or topological charge), which is a measure of the strength of the topological defect and is defined as the number of times that the order parameter turns of an angle of $2\pi$ while moving along a close contour surrounding the defect core. Hence possible values of defect strengths critically depend upon the nature of the order parameter: indeed polar systems only admit topological defects with integer winding numbers (Fig. 2b), while nematic systems offer a wider scenario; in fact by virtue of the head-tail symmetry, the headless nematic director can give rise to disclination patterns that also allows for half-integer winding numbers (Fig. 2a). Fig. 3 shows, for example, two defects of charge $\pm 1$ in an active contractile polar system: their mutual attracting interaction, due to elastic deformations, couples to the hydrodynamics generating a backflow Thot2002 ; giomi2013 that moves the two defects closer and leads to their annihilation. Fig. 3 also shows how defects act as a source of vorticity with the velocity field tilted with respect to polarization. On the contrary, if the system is extensile, activity drives defects of opposite topological charge apart and suppresses pair annihilation giomi2013 ; sanchez2012 . In simulations the correct position of a topological defect can be tracked either by looking at the polarization (or director for nematics) field profile or, only for nematics, by locating the regions where the scalar order parameter of the tensor field drops down. In the latter case, a further method, based on computing the degree of biaxiality around the defect core is briefly discussed in Appendix A. Indeed, regions close to the defect core display biaxiality PhysRevLett.59.2582 . Note that, although defects appearing in the active fluid of Fig. 3 are points, other structures are possible.

Defects are said to be topologically stable if a non uniform configuration of the order parameter cannot be reduced to a uniform state by a continuous transformation. A general criterion to establish whether a defect is topologically stable or not, is to look at the dimension $n$ of the order parameter. In a $d$-dimensional space, the condition that all the $n$ components of the order parameter must vanish at the defect core defines a “surface” of dimension $d-n$. Hence defects exist if $n\leq d$. In Fig. 3, for example, we have a two-dimensional system ($d=2$) with an order parameter (the polarization $\vec{P}$) having two components ($n=2$), and the defects allowed are points (or vortices). However, point defects can be unstable in quasi-$2d$ systems, i.e. when the order parameter fully lives in the three-dimensional space, as in such case one would have $n>d$: indeed the vector field in proximity of a vortex is always capable to escape out of the plane aligning with itself, thus removing the defect. In three-dimensional systems ($d=3$) one may have either point defects (if $n=3$) or lines (if $n=2$).

2.2 Free energy

In this Section we will shortly review the free-energy expressions generally used to describe polar and nematic suspensions and often employed in studying active fluids, built from the order parameters previously discussed.

Bulk properties and order-disorder phase transitions can be derived by a free energy functional with terms respecting the symmetries of the disordered phase, in the spirit of Landau approach. Free-energy $F$ will only contain scalar terms invariant under space rotations, proportional to the order parameters and their powers. For a vectorial order parameter, scalar objects of the form $\vec{P}^{2m}$ can be considered, with $m$ positive integer, usually arresting the expansion to the fourth order. For the nematic order parameter scalar quantities are of the form $\operatorname{Tr}(\vec{\underline{\underline{Q}}}^{m})$; note that there is no impediment here to odd power terms, by virtue of the invariance of $\vec{\underline{\underline{Q}}}$ under inversions, but no linear term will appear in the expansion since $\operatorname{Tr}\vec{\underline{\underline{Q}}}$ is identically null by definition. The presence of a third order term will lead to a first order nematic-isotropic transition through the establishment of metastable regions in the phase diagram chaikin1995 . Table 1 summarizes the bulk contributions to free energy for both polar and nematic systems; note that the uniaxial free energy can be derived from the biaxial case by writing the $\vec{\underline{\underline{Q}}}$ tensor through Eq. (4).

In order to take into account the energetic cost due to continuous deformations of the order parameters, elastic terms are also introduced in the free energy functional. In both polar and nematic systems three different kinds of deformations can be identified: splay, twist and bending, gauged to the theory through (in general) different elastic constants $\kappa_{1},\kappa_{2},\kappa_{3}$. While splay is related to the formation of fan-out patterns of the director and polarization field, bending generates rounded circular patterns. Instabilities associated to such deformations underlie the establishment of defects of different strength. Twist is forbidden in pure bidimensional systems, since this kind of deformation implies the director to coil around an axis, normal to the director itself. Table 1 also provides a picture of the energetic cost due to different kinds of deformations in terms of $\vec{P}$ and $\vec{n}$, respectively for polar systems and uniaxial nematics, under the assumption of uniform ordering ($S=cost$). The most general case is provided by the elastic contributions in biaxial nematics and still applies to the uniaxial case with $S=S(\vec{r})$. In order to exploit which terms are related to which deformations, one should expand the $\underline{\underline{\vec{Q}}}$ tensor into the elastic biaxial free energy in terms of the director through Eq. (4); doing so and grouping splay, twist and bend contributions one finds, after some algebric effort, that

[TABLE]

given that the Frank constants $\kappa_{i}$ fullfill the condition $\kappa_{3}\geqslant\kappa_{1}\geqslant\kappa_{2}$ to guarantee the positivity of $L_{i}$ schiele1983 . In many practical situations it is convenient to adopt the single constant approximation, consisting in setting all elastic constants equal to the same value, leading to a much simpler form for the elastic free energy chaikin1995 .

2.3 Active Forces

So far we reviewed the well known theoretical description for liquid crystals and fluids with anisotropic order parameter. We will see now how the active behavior of the constituents of the fluid can be expressed into the theoretical framework. The most direct way to develop the equations of motion for active systems at continuum level is by explicitly coarse-graining more detailed particle-based models ramaswamy2003 ; marchetti2013 . Therefore, before starting the theoretical description, we spend few words in describing the swimming mechanism of some microorganisms.

In general, the propulsive motion of active agents dispersed in a fluid creates a circulating flow pattern around each swimmer. The specific swimming mechanism of bacteria, for example, causes fluid to be expelled both forwards and backwards along the fore-aft axis, and drawn inwards radially towards this axis, creating an extensile flow pattern (Fig. 4). In some cytoskeleton extracts (such as the actomyosin protein complex), motor proteins can pull the filaments amongst themselves, causing them to contract lengthwise and giving rise to a contractile flow opposite to that of the previous example (Fig. 4)111A more detailed description of the hydrodynamics of swimmers is given in marchetti2013 ; ramaswamy2010 ; yeomans2017 .. Typically, activity creates a flow pattern that can be complicated in the near field, but whose far field is generically equivalent, at the lowest order, to the action of a force dipole pedley1992 and can be represented as such. By summing the contributions from each force dipole and coarse-graining simha2002 , it is possible to show that the stress exerted by the active particles on the fluid has the form

[TABLE]

where $\zeta$ is a phenomenological parameter that measures the activity strength, being negative for contractile systems and positive for extensile ones, while $\phi$ represents the concentration of the active material. Usually only terms linearly proportional to $\zeta$ are considered. In the case of polar active liquid crystals, the description can be carried out considering only the polarization field, re-expressing $\underline{\underline{\vec{Q}}}$ as a function of $\vec{P}$. The active stress in terms of the dynamical variable $\vec{P}(\vec{r},t)$ takes the form

[TABLE]

The expressions Eq. (6)) and Eq. (7), as we will see later, have been largely applied in the study of active fluids.

Many biological systems also display a local chirality Zhou2014 ; Naganathan2014 . Actin filaments, for example, are twisted in a right-handed direction Depue1965 so that myosin motors tend also to rotate them while pulling, creating a torque dipole. A concentrated solution of DNA has long been known to exhibit a cholesteric or blue-phase in different salt conditions livolant1986 ; livolant1991 . Such system can be made motile if interacting with DNA- or RNA-polymerases or with motor proteins. The effect of chirality, more than being taken into account by a suitable cholesteric term in the free energy222Chirality can be modelled wright1989 introducing a suitable term $(\nabla\times\underline{\underline{\mathbf{Q}}}+2q_{0}\underline{\underline{\mathbf{Q}}})^{2}$ in the free energy. This contribution favors the formation of an helix in the director pattern, whose pitch $p_{0}=2\pi/q_{0}$., can be incorporated in the description adding to the active stress extra terms, providing a source of angular momentum. For instance, if the active particles act on the surrounding fluid with a net torque monopole, a coarse-graining procedure furthauer2012 shows that a suitable choice for the nematic chiral stress tensor is given by $\zeta_{2}\epsilon_{\alpha\mu}Q_{\mu\beta}$ maitra2019 , where $\epsilon_{\alpha\mu}$ is the second order Levi-Civita tensor. Analogously, if the net torque is null but torque dipoles do not vanish, the corresponing stress tensor is given by $\zeta^{\prime}_{2}\epsilon_{\alpha\beta\mu}\partial_{\nu}\phi(P_{\mu}P_{\nu})$ Tjhung2017 , with $\epsilon_{\alpha\beta\mu}$ the third order Levi-Civita tensor. The sign of the second activity parameter $\zeta_{2}$ or $\zeta^{\prime}_{2}$ determines whether the stress generates a flux parallel ($\zeta_{2},\zeta^{\prime}_{2}>0$) or antiparallel ($\zeta_{2},\zeta^{\prime}_{2}<0$) with respect to the helicity of the twisting deformation. These terms drive the system out of equilibrium by injecting energy into it, and, as those of Eq. (6) and Eq. (7), cannot be derived from a free energy functional. In this approach the active stress tensor enters the hydrodynamic equations governing the motion of the self-propelled particles suspension, as discussed in Section 2.5. These are constructed from general principles, by assuming that an active gel is described by (a) “conserved” variables, which are the fluctuations of the local concentration of suspended particles and the total (solute plus solvent) momentum density, and (b) “broken-symmetry” variables, which, in the nematic phase, is the deviation of the director field from the ground state.

A more general way to construct the equations of motion at a coarse-grained level, is to generalize the forces-and-fluxes approach Degroot2013 to active systems kruse2004 . Considering for example an active gel characterized by polarization $\vec{P}$ and velocity $\vec{v}$, or equivalently by the strain rate tensor $u_{\alpha\beta}=(\partial_{\alpha}v_{\beta}+\partial_{\beta}v_{\alpha})/2$, the generalized hydrodynamic equations can be derived using Onsager relations, thus expanding fluxes $\partial_{t}\vec{P}$ and the stress tensor in terms of their conjugate forces $-\delta{F}/\delta\vec{P}$ and $u_{\alpha\beta}$ respectively, with $F$ polarization free energy. Active dynamics is obtained holding the system out of equilibrium by introducing a further pair of conjugate variables, namely the chemical potential difference between ATP and hydrolysis products and the rate of ATP consumption kruse2004 . This approach can be further generalized 1742-5468-2017-5-054002 including thermal fluctuations, recasting the forces-and-fluxes approach in the language of coupled generalized Langevin equations mazenko2006nonequilibrium .

Finally, we mention a more phenomenological model used to show self organization and scale selection for the flow pattern in active matter. This approach is inspired by the use of the Brazovskii model braz1975 ; Gonnella1998 for describing system with periodic order parameter and by dynamical approaches in studies regarding the role of hydrodynamics Swift1977 in the onset of convection Rayleigh-Bénard instability chandrasekhar1981hydrodynamic . Higher order derivatives of the velocity gradients are considered in the stress tensor in addition to the usual dissipative terms:

[TABLE]

If $\Gamma_{2}$ is chosen to be negative, this corresponds to the injection of energy in a definite range of wavelengths, while $\Gamma_{4}>0$ corresponds to hyperviscosity flow damping. This is obtained by truncating a long-wavelength expansion of the stress tensor bellout2013incompressible . The resulting generalized Navier-Stokes equations have been proven to capture experimentally observed bulk vortex dynamics of bacterial suspensions and some rheological properties of active matter dunkel2013_2 ; dunkel2013 ; wensink2012 .

2.4 Fluid mixtures with an active component

The active stress expressions of Eqs. (6) and (7) depend on the concentration of the active material. This quantity in turn can be a dynamical field if one would like to take into account a inhomogeneous presence of the active material in the solution. At level of particle description, different kinds of models for mixtures of self propelled and passive units have been considered. For example, Brownian-like simulations mccandlish2012 ; grosberg2015 ; stenhammar2015 focused on the role of activity in separating the two components of the mixtures. In a continuum description, binary fluids with an active component have been studied in tjhung2011 ; tjhung2012 ; nature2015 ; blow2014 showing that the active part may cause instabilities on an active-passive interface. Here we only introduce, as an example among the different models that can be used to describe fluid mixtures with an active component, the free-energy for a binary mixture where the active component is a polar gel tjhung2011 . It is given by

[TABLE]

The first term, multiplied by the phenomenological constant $a>0$, describes the bulk properties of the fluid; it is chosen in order to create two free-energy minima, one ($\phi=0$) corresponding to the passive material and the other one ($\phi=\phi_{0}$) corresponding to the active phase. The second one determines the interfacial tension between the passive and active phase, with $k$ positive constant. The third and the fourth terms control the bulk properties of the polar liquid crystal. Here $\alpha$ is a positive constant and $\phi_{\textrm{cr}}=\phi_{0}/2$ is the critical concentration for the transition from isotropic $(|\vec{P}|=0)$ to polar $(|\vec{P}|>0)$ states. The choice of $\phi_{\textrm{cr}}$ is made to break the symmetry between the two phases and to confine the polarization field in the active phase $\phi>\phi_{\textrm{cr}}$. The term proportional to $(\nabla{\vec{P}})^{2}$ describes the energetic cost due to elastic deformations in the liquid crystalline phase (see Table 1) in the single elastic constant approximation. Finally, the last term is a dynamic anchorage energy and takes into account the orientation of the polarization at the interface between the two phases. If $\beta\neq 0$, $\vec{P}$ preferentially points perpendicularly to the interface (normal anchoring): towards the passive (active) phase if $\beta>0$ ($\beta<0$). This choice for the anchoring is suggested by experimental observations. For instance, bacterial orientation at water-oil interfaces results from a relatively hydrophobic portion of each cell being rejected from the aqueous phase of the system Marshall1973 .

Such model can be also extended to study active nematic gels, by using the nematic tensor in place of the polarization field blow2014 ; giomi2013 ; giomi2014 ; doostmohammadi2018 . In this case the coefficients of the expansion of $Tr(\vec{\underline{\underline{Q}}}^{n})$ in bulk free energy (see Table 1) would depend on the scalar field $\phi$ and the elasticity, again written in the single elastic constant approximation, would include a term of the form $L\partial_{\alpha}\phi Q_{\alpha\beta}\partial_{\beta}\phi$ (with $L$ constant) to guarantee a perpendicular anchoring of the liquid crystal at the interface.

We finally mention a recent generalization of such models where emulsification of the active component is favored by the presence of surfactant added to the mixture bone2017 . This is done by allowing negative values of the binary fluid elastic constant $k$ and by including a term of the form $\frac{c}{2}(\nabla^{2}\phi)^{2}$ (with $c$ positive constant) to guarantee the stability of the free-energy.

A different continuum model, specifically introduced to study the motility induced phase separation (MIPS) without direct appeal to orientational order parameters $\vec{P}$ or $\underline{\underline{\vec{Q}}}$, but only to the scalar concentration field $\phi$, is the so called Active-model H tiribocchi2015 . In the old classification by Hohenberg and Halperin hohenberg1977 , the passive model H considers a diffusing, conserved, phase separating order parameter $\phi$ coupled to an isothermal and incompressible fluid flow through the advection-diffusion equation that will be introduced in Section 2.5. The chemical potential that enters the dynamic equation of the passive model H is given by

[TABLE]

with $a$, $b$, $k$ constants appearing in the Landau free energy for binary mixtures bray1994 (with $a$ negative in order to have phase separation between the two fluid components and $b$ and $k$ positive for stability). The same terms appear in Eq. (2.4) without the polarization contributions. The active model is then constructed by adding a leading order time-reversal breaking active term of the form $\mu_{a}=\lambda_{a}\left(\nabla\phi\right)^{2}$ (with $\lambda_{a}$ constant), not stemming from the free energy functional tiribocchi2015 . The deviatoric stress $\underline{\underline{\sigma}}$, that enters in the NS equations for the fluid flow, is, in $d$ dimensions,

[TABLE]

and can be obtained from the free energy, according to the formula reported in the second row of Table 2, only if $\hat{\zeta}=k$. If $\hat{\zeta}\neq k$ this is not true anymore and Eq. (11) is the sole leading-order contribution to the deviatoric stress for scalar active matter. Again here, $\zeta<0$ describes contractile systems while $\zeta>0$ the extensile ones. While $\mu_{a}$ has been found to create a jump in the thermodynamic pressure across interfaces and to alter the static phase diagram wittkowski2014 , the active stress $\underline{\underline{\sigma}}^{a}$ creates a negative interfacial tension in contractile systems that arrests the coarsening tiribocchi2015 .

2.5 Hydrodynamic equations

We can now introduce the hydrodynamic equations for active liquid crystals. Evolution equations for mass density $\rho(\vec{r},t)$ and velocity $\vec{v}(\vec{r},t)$ are given by

[TABLE]

with the energy balance equation generally neglected in this context. Eq. (12) is the continuity equation for mass density. In most of active matter systems Mach numbers $Ma$, defined as the ratio of the stream velocity and the speed of sound, is small; in such limit, this equation reduces to the solenoidal condition for the velocity field

[TABLE]

so that the fluid in this regime can be assumed at all practical effects as incompressible. Eq. (13) is the Navier-Stokes equation, where $p$ is the ideal fluid pressure and $\underline{\underline{\sigma}}$ is the stress tensor beris1994 that can be split into the equilibrium/passive and non-equilibrium/active contributions:

[TABLE]

The passive part is, in turn, the sum of three terms:

[TABLE]

The first term is the viscous stress, written as $\sigma_{\alpha\beta}^{viscous}=\eta(\partial_{\alpha}v_{\beta}+\partial_{\beta}v_{\alpha})$, where $\eta$ is the shear viscosity333In the compressible case, the viscous stress tensor also includes a term proportional to the divergence of the velocity, such that:

$\sigma_{\alpha\beta}^{viscous}=\eta(\partial_{\alpha}v_{\beta}+\partial_{\beta}v_{\alpha})+\left(\tilde{\zeta}-\frac{2\eta}{d}\right)\partial_{\gamma}v_{\gamma}\delta_{\alpha\beta},$

(17)

where we denoted the bulk viscosity with $\tilde{\zeta}$. . An explicit form for the elastic and interface stress is reported for the polar and nematic cases in Table 2.

The order parameter $\mathbf{\Psi}$ of the active liquid crystal (that is $\underline{\underline{\vec{Q}}}$ for nematics and $\vec{P}$ for polar systems) evolves according to

[TABLE]

known as Beris-Edwards equation, within the theory of liquid crystal hydrodynamics described through the $\vec{\underline{\underline{Q}}}$-tensor. The term $\vec{S}$ accounts for the response of the orientational order to the extensional and rotational components of the velocity gradient and is reported for the polar stark2003 ; marchetti2006 and nematic degennes1993 case in the fourth row of Table 2. The molecular field $\mathbf{\Xi}$ governs the relaxation of the orientational order to equilibrium, and is multiplied by a collective rotational-diffusion constant $\Gamma$. Its expressions are given in the third row of Table 2. The left-hand side of Eq. (18) is commonly addressed as material derivative of the order parameter $\mathbf{\Psi}$, and can be formally derived making use of Liouville equations. In fact one can write $D_{t}\mathbf{\Psi}=\partial_{t}{\mathbf{\Psi}}+\{\mathbf{\Psi},\mathcal{H}\}$, where $\{...\}$ are the Poisson brackets and the Hamiltonian is $\mathcal{H}=F+\frac{1}{2}\int\rho\vec{v}^{2}$.

A more phenomenological procedure to derive the material derivative explicitly is based on the fact that order parameters can be advected by the fluid. Here we outline the procedure referring only to the polarisation field. We first note that the relative position $\vec{\tilde{r}}$ of two close points in the fluid evolves according to the following equation:

[TABLE]

where $\underline{\underline{D}}$ and $\underline{\underline{\Omega}}$ have been defined in Table 2. The first two contributions are the usual lagrangian derivative terms, while the third and fourth ones account respectively for rigid rotations and deformations of the fluid element. Thus the material derivative for the polarisation field will include the first three terms since a vector advected by the flow is capable to follow any rigid motion; for what concerns the last term in Eq. (19), this cannot enter directly into the material derivative of a vector field, but it must be weighted through an alignment parameter $\xi$, ruling the dynamical behavior of the vector field under enlargement and/or tightening of flow tubes. This allows us to obtain the material derivative for the polarization field simply substituting $\vec{P}$ in place of $\vec{\tilde{r}}$.

Finally the time evolution of the concentration field $\phi(\vec{r},t)$ of the active material is governed by an advection-diffusion equation

[TABLE]

where $M$ is the mobility and $\delta F/\delta\phi$ is the chemical potential. A more generalized form of the material derivative has been used to model self advective phenomena, for example, actin polymerization in motile eukaryotic cells tjhung2012 , by substituting $\nabla\cdot(\phi\vec{v})\rightarrow\nabla\cdot(\phi\vec{v}+w\vec{P})$, where $w$ is a constant related to the velocity of actin polymerization.

3 Lattice Boltzmann Method

A certain number of approaches are feasible when dealing with the description of fluid systems; each of them can be classified according to the level of spatial approximation. A molecular approach would hardly access the time and space scales relevant for a complete hydrodynamic description of the systems here considered. At a mesoscopic level, kinetic theory furnishes a description of irreversible and non-equilibrium thermodynamic phenomena in terms of a set of distribution functions encoding all necessary informations related to space positions and velocities of particles. Continuum equations give a description of irreversible phenomena by using macroscopic variables slowly varying in time and space. This last approach has the not-negligible advantage that one has to deal with a few fields. On the other hand, when considering continuous equations, one has to face some technical issues arising from the stability of numerical implementation and discretization schemes LeVeque2007 . Moreover, many numerical methods aimed at solving the continuous equations, exhibit criticalities in the amount of computational resources, mostly in terms of processing times and memory requirement, or in the implementation of boundary conditions in complex geometries. To avoid these issues lattice-gas-automaton (LGA) models were first developed starting from the pioneering work of Frisch et al. Frisch1986 . This kinematic approach to hydrodynamics is based on the description of the dynamics of a number of particles moving on a suitable lattice. An exclusion principle is imposed to restrict the number of particles with a given velocity at a certain lattice point to be [math] or $1$. This latter feature allows for a description of the local particle equilibrium through the Fermi-Dirac statistics alexander1992 . Despite LGA proved to be very efficient in simulating the Navier-Stokes equation from a computational point of view and in managing boundary conditions, LGA simulations are intrinsically noisy due to large fluctuations of local density. Moreover, they suffer from non-Galilean invariance, due to density dependence of the convection coefficient and from an unphysical velocity dependence of the pressure, arising directly from the discretization procedure Chen1992 .

Lattice Boltzmann methods were then developed to overcome these difficulties higuera1989 . Particles in the LGA approach are formally substituted by a discretized set of distribution functions, so that hydrodynamic variables are indeed expressed at each lattice point in terms of such distribution functions. Despite the fact that lattice Boltzmann is a mesoscopic numerical method, it has a number of advantages that resulted in a broad usage in many branches of hydrodynamics. Firstly LB algorithms are appreciably stable and they are characterized by their simplicity in the treatment of boundary conditions. Not to be neglected is the fact that LB algorithms are particularly suitable to parallel approach.

In the following of this Section we will first provide a simple overview of the method, without getting too technical, in order to convey to the reader the purpose of this approach. In Section 3.1 we will first introduce LBM for a simple fluid, while Section 3.2 will be devoted to recover the continuum hydrodynamic equations, already presented in Section 2.5. In Section 3.3 we will describe some routes to adapt LBM to the case of complex fluids by introducing either a forcing term or properly fixing the second moment of the equilibrium distribution functions. In Section 3.4 different approaches to deal with the advection-relaxation equations for order parameters coupled to the momentum equations will be examined. In Section 3.5 we will focus on some algorithms that have been recently used in the numerical investigation of active matter. Finally, in Section 3.6 we will focus on the computational performance and stability of the method.

3.1 General features of lattice Boltzmann method

The lattice Boltzmann approach to hydrodynamics is based on a phase-space discretized form of the Boltzmann equation Benzi1992 ; Rothman1997 ; Chopard1998 ; Wolf-Gladrow2000 ; succi2001 for the distribution function $f(\vec{r},\vec{\xi},t)$, describing the fraction of fluid mass at position $\vec{r}$ moving with velocity $\vec{\xi}$ at time $t$. Since space and velocities are discretized, the algorithm is expressed in terms of a set of discretized distribution functions $\{f_{i}(\vec{r}_{\alpha},t)\}$, defined on each lattice site $\vec{r}_{\alpha}$ and related to a discrete set of $N$ lattice speeds $\{\vec{\xi}_{i}\}$, labelled with an index $i$ that varies from $1$ to $N$ (see Fig. 5). In the case of the collide and stream version of the algorithm, the evolution equation for the distribution functions has the form

[TABLE]

where $\mathcal{C}(\{f_{i}\},t)$ is the collisional operator that drives the system towards equilibrium, represented by a set of equilibrium distribution functions, and depends on the distribution functions; its explicit form will depend upon the particular implementation of the method. Eq. (21) describes how fluid particles collide in the lattice nodes and move afterward along the lattice links in the time step $\Delta t$ towards neighboring sites at distance $\Delta x=\xi_{i}\Delta t$. This latter relationship is no more considered in finite difference lattice Boltzmann models (FDLBM) Cao1997 ; Mei1998 ; lee2001 ; Sofonea2003 ; Sofonea2004 ; Cristea2006 ; Cristea2010 . In this kind of models the discrete velocity set can be chosen with more freedom, making possible to use non uniform grids, selecting lattice velocities independently from the lattice structure444When dealing with FDLBM it is useful to introduce more than only one set of distribution functions $\{f_{ki}\}$, where the extra index $k$ labels different sets of discrete velocities $\{\vec{\xi}_{ki}\}$, with index $i$ still denoting the streaming direction. The evolution equation for distribution functions for the FDLBM reads:

$\partial_{t}f_{ki}+(\vec{\xi}_{ki}\cdot\nabla)f_{ki}=\mathcal{C}(\{f_{ki}\},t).$

(22)

Here differential operators must be discretized: Runge-Kutta or midpoint schemes can be used to compute the time derivative while there are several possibilities to compute the advective term on the left-hand side of the previous equation. For the reader interested in details of the implementation we suggest to refer to Sofonea2003 ; watari2003 .. This result is found to be extremely useful when it is necessary to release the constraint of having a constant temperature in the system watari2003 ; gonnella2007 . Moreover it might be also helpful in the case of LB models for multicomponent systems where the components have different masses and this would result in having different lattice speeds, one for each fluid species. Beside the wider range of applicability of the FDLBM with respect to the LBM, the latter furnishes a simple and efficient way to solve hydrodynamic equations; in addition we are not aware of any implementation of the FDLBM algorithm developed to study active matter; for this reasons we will avoid any further discussion on this variant of LBM.

In the case of a simple fluid, in absence of any external force, assuming the BGK approximation with a single relaxation time bhatnagar1954 , one writes

[TABLE]

where $f_{i}^{eq}$ are the equilibrium distribution functions and $\tau$ is the relaxation time, connected to the viscosity of the fluid, as it will be seen. The mass and momentum density are defined as

[TABLE]

where summations are performed over all discretized directions at each lattice point. By assuming both mass and momentum density to be conserved in each collision, it is found that conditions in Eq. (24), (25) must hold also for the equilibrium distribution functions:

[TABLE]

Moreover, it is necessary to introduce further constraints on the second moment of the equilibrium distribution functions to recover continuum equations, as it will become more evident in the following. Further constraints on higher order moments may become necessary to simulate more complex systems: for instance full compressible flows or supersonic adaptation of the algorithm may require the specification of moments up to the third, while for a complete hydrodynamic description in which heat transfer is also taken into account, even the fourth moment needs to be specified watari2003 . Active matter systems such as bacterial and microtubules suspensions reasonably fulfil the incompressible condition, so that in the following we will only impose constraints up to second order moments.

Another peculiar fact is that viscosity explicitly depends on the choice of a particular lattice. Due to the fact that sufficient lattice symmetry is required to recover the correct Navier-Stokes equation in the continuum limit Frisch1986 , not all the possible lattice structures can be adopted. By denoting the space dimension by $d$ and the number of lattice speeds by $Q$, Table 3 shows the velocities $\{\vec{\xi}_{i}\}$ and the corresponding weights in the equilibrium distribution functions (see next Section) for the most frequent choices. Here the quantity $c=\Delta x/\Delta t$, connected to the speed of sound of the algorithm, has been introduced as the ratio between the lattice spacing $\Delta x$ and the time step $\Delta t$. Figure 5 explicitly illustrates the lattice structures in the two-dimensional case.

3.2 Lattice Boltzmann for a simple fluid

In this Section we will present a basic lattice Boltzmann algorithm to solve the hydrodynamic equations (12) and (13) for a simple fluid; in this case the term on the right hand side of the Navier-Stokes equation (13) reduces to the pressure gradient plus the mere viscous contribution $\partial_{\beta}\sigma_{\alpha\beta}^{viscous}$, if no external force is acting on the fluid.

Conditions (26) and (27) can be satisfied by expanding the equilibrium distribution functions up to the second order in the fluid velocity $\vec{v}$ Chen1992 :

[TABLE]

where index $s=|\vec{\xi}_{i}|^{2}/c^{2}$ relates the $i$-th distribution function to the square module of the corresponding lattice velocity, and the greek index denotes the Cartesian component. This expansion is valid as far as the Mach number $Ma=v/c_{s}$ is kept small, $c_{s}$ being the speed of sound, whose explicit expression in turn depends upon the lattice discretization shan2006 . The present assumption has the important consequence that LB models based on the previous expansion of the equilibrium distribution functions have great difficulty in simulating compressible Euler flows, that usually take place at high Mach numbers. This issue arises in standard LB approaches because of the appearence of third order nonlinear deviations from the Navier-Stokes equation guangwu1998 . Qian and Orzsag demonstrated in qian1993 that such nonlinear deviations grow together with $Ma^{2}$, so that they can be neglected in the low Mach number regime but become important in the compressible limit555In order to overcome the limit posed by the low Mach number regime, many variations of the standard LBM have been developed. Alexander et al. proposed a model where the high Mach number regime could be achieved by decreasing the speed of sound alexander1992 ; discrete-velocity models nadiga1995 ; sun1998 were later introduced allowing for simulation of the compressible Euler equation in a wider range of Mach numbers. Other implementations are based on a Taylor expansion of the equilibrium distributions up to higher orders together with suitable constraints on the third and fourth moments shan2006 ; kataoka2004 ; kataoka2004_2 . . For such reasons it is necessary to ensure that velocities never exceed a critical threshold that can be reasonably chosen such that $Ma\lesssim 0.3$ qian1993 .

Besides constraints expressed by Eq. (26) and (27), an additional condition on the second moment of the equilibrium distribution functions is imposed so that

[TABLE]

This is a necessary condition to recover the Navier-Stokes equation in the continuum limit. By substituting the expansion in Eq. (28) in constraints introduced in Eq. (26), (27) and (29), a suitable choice for the expansion coefficients is found to be

[TABLE]

where for the sake of clarity we have explicitly chosen a $d2Q9$ lattice geometry. Requiring suitable isotropy conditions and Galilean invariance Chen2008 , it is even possible to show analytically Qian1992 , that the equilibrium distribution functions can be written in a more general way as

[TABLE]

where the weights $\omega_{i}$ are given in Table 3. In Appendix B it will be shown that the algorithm here presented correctly reproduces Eqs. (12) and (13).

We add for completeness that it is also possible to adopt a discretization in velocity space based on the quadrature of a Hermite polynomial expansion of the Maxwell-Boltzmann distribution shan2006 . One then gets a lattice Boltzmann equation that allows us to exactly recover a finite number of leading order moments of the equilibrium distribution functions. In this case the quantity $c$ is fixed and given by $c=2$ for the geometry $d2Q7$ and by $c=\sqrt{3}$ for the other geometries in Table 3. For a detailed discussion the interested reader may refer to Ref. shan2006 . Finally, we mention that it would be possible to introduce small thermal fluctuations into the algorithm, in a controlled way, by means of a stochastic collision operator. The fluctuation-dissipation theorem can then be satisfied by requiring consistency with fluctuating hydrodynamics Adhikari2005 . Since to the best of our knowledge there are no LB models for active systems including thermal noise, we do not give further details referring the interested reader to the Ref. Dunweg2009 .

3.3 LBM beyond simple fluids

So far we have implemented a lattice Boltzmann method for a simple fluid in absence of any forcing term, with only viscous contribution to the stress tensor. On the other hand when dealing with more complex systems, such as multicomponents or multiphase fluids, the stress tensor may include further contributions (such as elastic and interfacial ones, see Table 2) which have a non-trivial dependence on order parameters and their derivatives. In this Section we will show two different strategies adopted to numerically implement such terms. Briefly, while in the first one they are included in an extra term, appearing in the second moment of the equilibrium distribution functions, in the second one they enter through an external forcing added to the collision operator in the lattice Boltzmann equation.

3.3.1 First method

To implement a general symmetric stress tensor contribution in the lattice Boltzmann scheme previously introduced, we again impose the constraints of Eq. (26) and Eq. (27) on the zeroth and on the first moment of the equilibrium distribution functions, while constraint on the second moment previously given in Eq. (29) is modified according to the following relation

[TABLE]

Here $\sigma_{\alpha\beta}$ stands for the total stress tensor including preassure contributions, but deprived of viscous ones. Note that, due to the symmetry of the left hand side of Eq. (35), this algorithm can be applied to models that involve only symmetric contributions to the stress tensor. For instance, this method is suitable to study binary mixtures, as the stress tensor associated to the concentration contribution is indeed symmetric, but not liquid crystals, as the antisymmetric part of the relative stress tensor does not vanish (see Table 2). This latter case will be discussed in the following Sections. To satisfy Eq. (35) denniston2001 ; orlandini2008 , the equilibrium distribution functions can be expanded as follows

[TABLE]

where an extra term, quadratic in lattice velocities, has been added with respect to the case of a simple fluid (see Eq. (28)), to include a general stress tensor in the model. As for a simple fluid, the coefficients of the expansion can be calculated by imposing constraints of Eq. (26), Eq. (27) and Eq. (35). For a $d2Q9$ geometry a suitable choice is given by

[TABLE]

where we denoted by $\sigma_{\alpha\beta}^{0}$ the traceless part of $\sigma_{\alpha\beta}$.

One can now proceed to recover the Navier-Stokes equation by using a Chapman-Enskog expansion666The second moment constraint on the equilibrium distribution functions is not necessary for the derivation of the continuity equation. Hence the procedure to recover this equation is not affected by the modifications introduced in the new version of the algorithm, with respect to the case of a simple fluid.. Assuming that the fluid is flowing at small Mach numbers, so to ignore third-order terms in the fluid velocity, and taking the first moment of Eq. (77), one gets

[TABLE]

which is the Navier-Stokes equation at first order in Knudsen number. To recover the Navier-Stokes equation at second order, we start from Eq. (91), where we need to evaluate the second moment of $f_{i}^{(1)}$

[TABLE]

The first time derivative in square brackets is negligible at the leading order, while

[TABLE]

that shows, together with Eq. (38), that this term gives a null contribution. Finally, using Eq. (93) we get the same result of Eq. (95) which allows one to restore the Navier-Stokes equation.

3.3.2 Second method

An alternative route to the solution of the LB equation (21) relies on the use of a pure forcing method li2007 ; tiribocchi2009 . In this case the total stress tensor enters the model via a forcing term $\mathcal{F}_{i}$, without any additional constraint on the second moment of the equilibrium distribution functions, with condition given in Eq. (29). The collision term $\mathcal{C}_{f_{i}}$ assumes the simple form of the BGK approximation supplemented by a forcing term

[TABLE]

where the equilibrium distribution functions $f_{i}^{eq}$ are again expressed as a second-order expansion in the velocity $\vec{v}$ of the Maxwell-Boltzmann distribution Qian1992 . The fluid momentum is now given by the average between the pre- and post-collisional values of the velocity $\vec{v}$, as usually done when using a forcing term Shan1994 ; Buick2000

[TABLE]

where $F_{\alpha}$ is the cartesian component of the force density acting on the fluid. The choice of the equilibrium distribution functions and their constraints is kept as in Section 3.2, with coefficients given by Eqs. (30)-(33) for a $d2Q9$ lattice. The term $\mathcal{F}_{i}$ can be written as an expansion at the second order in the lattice velocity vectors ladd2001 :

[TABLE]

where coefficients $A$, $B_{\alpha}$ and $C_{\alpha\beta}$ are functions of $F_{\alpha}$. In order to correctly reproduce hydrodynamic equations, the moments of the force term must fulfil the following relations

[TABLE]

which lead to guo2002

[TABLE]

To recover the continuity (12) and the Navier-Stokes (13) equations it suffices to require that

[TABLE]

From the Chapman-Enskog expansion (see Appendix C for the details of the calculation) it results that the fluid viscosity in Eq. (17) is $\eta=\frac{\rho\Delta tc_{s}^{2}}{3}(\tau-1/2)$. No extra contributions appear in the continuum equations (12) and (13), apart from a term of order $v^{3}$ which can be neglected if the Mach number is kept small.

Other approaches to the numerical solution of the LB equation introduce spurious terms which cannot always be kept under control. For a complete discussion the interested reader may refer to Ref. guo2002 . The one presented here has proved to be effective for simple fluids guo2002 , multicomponent gonnella2010 and multiphase fluid systems coclite2014 ; sofonea2018 even though, as far as we know, a full external forcing algorithm has not been applied to active systems yet. We add that boundary walls can be easily implemented as illustrated in the Appendix D.

3.4 Coupling with advection-diffusion equation

The aim of lattice Boltzmann methods goes far beyond the treatment of Navier-Stokes equation; indeed, it has proven to be a fundamental tool to solve general conservation equations ancona1994 . Moreover, beside many implementations devoted to hydrodynamics studies, such as the ones cited at the end of the previous Section, recently a LBM approach has also been used to solve Einstein equations for gravitational waves Ilseven2016 .

We devote this Section to report on two characteristic ways to solve the dynamics of order parameters coupled to hydrodynamics in a fluid system. Because of its relevance in the study of complex fluids we will focus on the treatment of the advection-diffusion equation (20) for a concentration field. The first possibility is to develop a full LBM approach in which the advection-diffusion equation is solved by introducing a new set of distribution functions $\{g_{i}(\vec{r},t)\}$ connected to the concentration field, beside the distribution functions $\{f_{i}(\vec{r},t)\}$ needed to solve the Navier-Stokes equation. Another route is to follow a hybrid approach where the advection-diffusion equation is solved via a standard finite difference algorithm while hydrodynamics is still solved through a LB algorithm.

Full LBM approach

To solve the hydrodynamic equations for a binary system through a full LB approach the introduction of a new set of distribution functions $\{g_{i}(\vec{r},t)\}$ is needed xu2003 ; xu2004 . The index $i$ again assigns each distribution function to a particular lattice direction indicated by the velocity vector $\vec{\xi}_{i}$. The concentration field $\phi(\vec{r},t)$ is thus defined as

[TABLE]

As in Eq. (21), distribution functions $g_{i}$ evolve according to the following equation

[TABLE]

where the BGK approximation for the collisional operator has been used. A new relaxation time $\tau_{\phi}$ has been introduced since the relaxation dynamics of the concentration field may consistently differ from that of the underlying fluid. In Eq. (48) we have also introduced the set of equilibrium distribution functions $\{g_{i}^{eq}(\vec{r},t)\}$ that fulfill the following relation

[TABLE]

This ensures that the concentration field is conserved during the evolution.

To recover the advection-diffusion equation in the continuum limit, it is necessary to impose the following constraints on the first and second moments of the equilibrium distribution functions

[TABLE]

Here the mobility parameter $\chi$ tunes the diffusion constant $M$ that appears on the right-hand side of the advection-diffusion equation, while $\mu$ is the chemical potential. A suitable choice of the distribution function which fulfills Eq. (49), Eq. (50) and Eq. (51) can be written as a power expansion up to the second order in the velocity

[TABLE]

where the coefficients of the expansion can be computed from Eqs. (37) through the formal substitution

[TABLE]

The continuum limit of the advection-diffusion equation can be performed through a Taylor expansion of the left-hand side of Eq. (48) and by using Eqs. (49)-(51) yeomans2009 . This leads to the following expression of the diffusion constant

[TABLE]

This algorithm can be generalized to describe the evolution of more complex order parameters, such as the nematic tensor $Q_{\alpha\beta}$, whose dynamics is governed by the Beris-Edwards equation of motion (Eq. (18)). Since $Q_{\alpha\beta}$ is a traceless symmetric tensor, in $d$ dimensions, at least $d(d+1)/2-1$ extra distribution functions $\{G_{i,\alpha\beta}(\vec{r},t)\}$ are needed, which are related to $Q_{\alpha\beta}$ through

[TABLE]

The rest of the algorithm can be thus developed as the one presented for the concentration field. In Section 3.5 we will go back to LBM for liquid crystal dynamics and we will present another algorithm that employs a predictor-corrector numerical scheme.

Hybrid LBM approach

An alternative approach to solve the Navier-Stokes equation and an advection-diffusion equation for an order parameter is based on a hybrid method, in which a standard LBM solves the former while a finite-difference scheme integrates the latter equation.

Let us consider, for instance, the evolution Eq. (20) of the concentration field $\phi(\vec{r},t)$. Space $\vec{r}$ and time $t$ can be discretized by defining a lattice step $\Delta x_{FD}$ and a time step $\Delta t_{FD}$ for which $\Delta x_{FD}=\Delta x_{LB}$ (namely the scalar field is defined on the nodes of the same lattice used for the LB scheme) and $\Delta t_{LB}=m\Delta t_{FD}$, with $m$ positive integer. At each time step the field $\phi$ evolves according to Eq. (20) and is updated in two partial steps.

Update of the convective term by means of an explicit Euler algorithm

[TABLE]

where all variables appearing at the right-hand side are computed at position $\vec{r}_{\alpha}$ and time $t$. Note that the velocity field $\vec{v}$ is obtained from the lattice Boltzmann equation. 2. 2.

Update of the diffusive part

[TABLE]

Note that one could use more elaborate methods to solve convection-diffusion equations. For instance, one can combine predictor-corrector schemes for the treatment of the advective term with a wealth of finite-difference schemes for the numerical solution of parabolic equations Evans2000 . Nevertheless one has to always keep in mind consistency between the order of accuracy of combined different numerical schemes used. However, the method here described, besides being relatively simple to implement, combines a good numerical stability with a reduced memory requirement with respect to the full LBM approach tiribocchi2009 , as it will be discussed in Section 3.6.

3.5 LBM for Active Fluids

As outlined in Section 2, many properties of active matter are captured by liquid crystal hydrodynamics. Here we describe a LB method that solves both the Navier-Stokes equation and the Beris-Edwards equation through a full LB approach, a method often employed to numerically investigate active matter denniston2001 ; Cates2099 .

As the liquid crystal stress tensor entering the Navier-Stokes equation is generally not symmetric, one could either (i) build an algorithm in which it is fully included through an external forcing term (as described in Section 3.3.2) or (ii) separate the symmetric part from the antysimmetric one, by including the former in the second moment of the equilibrium distribution functions and treating the latter as an external forcing term. Although the two procedures are equivalent, only the second approach, first introduced by Denniston et al. denniston2001 , has been developed so far.

In this method two sets of distribution functions, $\{f_{i}\}$ and $\{G_{i,\alpha\beta}\}$, are defined and are connected to the hydrodynamic variables (i.e. density, momentum and order parameter) through Eqs. (24), (25) and (55). Their evolution equations are solved by using a predictor-corrector-like scheme

[TABLE]

[TABLE]

where $f_{i}^{*}$ and $G_{i,\alpha\beta}^{*}$ are respectively first order approximation to $f_{i}^{*}(\vec{r}+\vec{\xi}_{i}\Delta t,t+\Delta t)$ and $G_{i,\alpha\beta}^{*}(\vec{r}+\vec{\xi}_{i}\Delta t,t+\Delta t)$ obtained by setting $f_{i}^{*}\equiv f_{i}$ and $G_{i,\alpha\beta}^{*}\equiv G_{i,\alpha\beta}$ in Eq. (58) and (59). The collisional terms are given by a combination of the usual collision operator in the BGK approximation plus a forcing term

[TABLE]

where $\tau_{f}$ and $\tau_{G}$ are two distinct relaxation times, and $p_{i}$ and $M_{i,\alpha\beta}$ are the two additional forcing terms.

In order to recover continuum equations one must impose constraints on the zeroth, first and second moments of the equilibrium distribution functions and on the forcing terms. The local conservation of mass and momentum is ensured by (26) and (27), while the second moment is given by Eq. (35), in which the stress tensor on the right hand side includes the sole symmetric part. The antisymmetric contribution $\sigma^{anti}_{\alpha\beta}$ is introduced through the forcing term $p_{i}$, which fulfills the following relations

[TABLE]

The remaining distribution functions $G_{i,\alpha\beta}$ obey the following equations

[TABLE]

while the forcing term $M_{i,\alpha\beta}$ satisfies

[TABLE]

We finally note that the predictor-corrector scheme has been found to improve numerical stability of the algorithm and to eliminate lattice viscosity effects (usually emerging from the Taylor expansion and appearing in the viscous term, in the algorithms discussed so far) to the second order in $\Delta t$. To show this, one can Taylor expand Eq. (58) to get

[TABLE]

The left-hand side is $\mathcal{O}(\Delta t)$ and coincides with the term in square brackets. One could then write at second order in $\Delta t$

[TABLE]

An analogous calculation for $G_{i,\alpha\beta}$ shows that

[TABLE]

thus recovering the proper lattice Boltzmann equations.

A hybrid version of the algorithm, widely employed in the study of active matter, solves the Navier-Stokes equation through a predictor-corrector Lattice-Boltzmann approach and the Beris-Edward equation by means of a standard finite-difference method henrich2010 ; Cates2099 .

Further models involving more than just one order parameter have been developed in recent years, such as the theory discussed in Section 2.4, in which the liquid crystal order parameter (the polarization field) is coupled to the concentration field of a binary fluid mixture. Again a hybrid approach, in which both equations of the concentration and of the polarization have been solved through finite difference methods, has been used bone2017 ; negro2018 .

3.6 Computational perspectives: stability, efficiency and parallelization

In the previous Sections we presented different LB algorithms for the treatment of the hydrodynamics of complex and active fluids. We will comment here on the stability of two different $d2Q9$ hybrid LB codes solving the equations of an active polar binary mixture (the hydrodynamics is solved by means of LB while the order parameter dynamics is integrated by a finite difference algorithm implementing first order upwind scheme and fourth order derivative accuracy), described by the free energy in Eq. (2.4), treating the symmetric part of the stress tensor with two different approaches. The first is a mixed approach, presented in the previous Section, where the symmetric part of the stress tensor enters in the definition of the second moment of the distribution functions (see Eq. (35)) while the anti-symmetric part is treated by means of the forcing term $p_{i}$ (see Eqs. (60) and (62)). In the second approach the total stress tensor is treated by means of the only forcing term. To compare the stability of the two algorithms we fixed the mesh spacing and the time resolution ($\Delta x=1$, $\Delta t=1$), and we let vary the relaxation time $\tau$ and the intensity of active doping $\zeta$ appearing in the active stress tensor (7). The results of the stability test in Fig. 6 show that the full-force approach is definitely more stable than the mixed one. In this latter case the code is found to be stable for $\tau>0.715$ in the passive limit ($\zeta=0$) while to simulate active systems ($\zeta>0$), the relaxation time must be accurately chosen to ensure code stability. In the full-force approach the code is found to be stable for $\tau>0.5$, almost independently of $\zeta$.

The rest of this Section is devoted to a brief discussion of some performance aspects, such as efficiency and parallelization of a LB code. LBM is computationally efficient if compared to other numerical schemes. The reason lies in the twofold discretization of the Boltzmann equation in the physical and velocity space. For instance, computational methods such as finite-difference (FD) and pseudo-spctral (PS) methods require high order of precision to ensure stability Evans2000 and to correctly compute non-linearities in the NS equation (13). This introduces non-local operations in the computational implementation that reduce the throughput of the algorithm. LBM, on the contrary, is intrinsically local, since the interaction between the nodes is usually more confined, according to the particular choice of the lattice, while non linearities of the NS equation is inherently reproduced at the level of the collision operator. For instance, while the number of floating point operations needed to integrate the hydrodynamics equations on a $d$-dimensional cubic grid is $\sim L^{d}$ for LBM, it is instead of order $\sim(\ln L)L^{d}$ for pseudo-spectral models succi1991 . Nevertheless LB algorithms are definitely much more memory consuming, since for each field to evolve, one needs a number of distribution functions equal to the number of lattice velocities. From this perspectives, the hybrid version of the code is somewhere in the middle between the two approaches, since it allows one to exploit both computational efficiency and simplicity typical of LB approaches and, at the same time, to keep the amount of memory to be allocated at runtime lower than that necessary for a full LB treatment.

LB algorithms are also suitable for parallelization. The reason still lies in the local character of LB, since at the base of the efficiency of any parallelization scheme is the compactness of the data that must be moved among the different devices that take part in the program execution. Parallelization approaches involving both CPUs, i.e. MPI or OpenMP, and GPUs, such as CUDA and OpenCL, or even both (CUDA aware MPI) can be used when dealing with LB kruger2016 . Most of them, such as OpenMP or GPU-based approaches, aim at rising the amount of floating operations per unit time, while a different technique consists in splitting the global computational domain in subdomains and assign each of them to a different computational unit (usually threads of one or more processors). This is usually done with MPI.

Fig. 7 shows the results of a strong scaling test performed on a hybrid code integrating the hydrodynamics of a polar binary mixture bone2017 ; negro2018 , implementing the full-force algorithm used for the stability analysis. This test consists in changing the amount of processors used to perform a certain task, while keeping fixed the size of the computational grid and measuring the speedup, namely the ratio of time spent to perform the operation with only one processor over the time taken when more processors are used. Simulations were performed both with $d2Q9$ (hollow dots) and $d3Q15$ (full dots) lattice structures on different computational infrastructures (Archer (red), Marconi (blue) and ReCas (green)). While for a few number of processors the scaling is approximately linear, thus close to the ideal linear behavior (black line), as the number of processors increases, it progressively deviates from the ideal scaling law. This is due to a number of issues that may depend both on the infrastructure characteristics (bandwidth, cache size, latency, etc.) and on the program implementation (bottlenecks, asynchrony among processor, etc.). Moreover, code scalability is found to be significantly better in three-dimensional grids than in their bidimensional version. This is because the fraction of time spent by the $3d$ code to perform parallel operations (sending and receiving data, reduction operations, synchronization, etc.) is consistently reduced compared with its $2d$ counterpart.

4 Spontaneous flow

Many remarkable phenomena in the physics of active fluids are related to the flow behavior induced by the presence of active forcing in the dynamical description of the system. The first effect that was studied, to which this section is devoted, is the occurrence of spontaneous flow in fluids with sufficiently strong activity. Numerical methods, and LBM in particular, have been essential in the study of this problem, both for supporting theoretical analysis of instabilities and for the understanding of non-linear regimes. Here we will put the accent on numerical studies, while a more complete review of analytical studies of instabilities in active fluids can be found in marchetti2013 .

The continuum theory described in Section 2.5 with Eqs. (6) or (7) was deduced in simha2002 where instabilities due to activity were first analyzed. The instability of ordered contractile and extensile gels to splay and bend perturbations respectively, giving rise to macroscopic flows giomi2008 ; edwards2009 , can be qualitatively explained as follows. Let us consider the case of contractile dipoles initially perfectly ordered, as in the sketch on the left of Fig. 8. In this situation the force dipoles balance each other and the net flow, obtained by the sum of those due to single dipoles as represented in Fig. 4, is null. However, if a small splay deformation is present (middle panel in Fig. 8), the density of contractile forces on the left is larger than that on the right, and a flow sets up. This flow causes further splay which destabilizes the system that starts to flow macroscopically. For extensile activity, under the same splay deformation, the initial flow (directed to the left in this case, see Fig. 4) would align the dipoles and no net macroscopic flow would appear. By a similar argument, it can be shown that extensile fluids get unstable to bend deformations.

The spontaneous flow instability for contractile systems was analyzed in voituriez2005 for the simple geometry of a bidimensional thin film confined on a one-dimensional substrate, with planar anchoring on the confining boundaries. For small thicknesses or small activity, boundary effects are prevailing and the gel remains in an unperturbed, static, homogeneously polarized state. Above a critical thickness or a critical activity, a polarization tilt appears and the system flows with a finite shear gradient. The study has been later extended to films where undulations of the free surface are also considered sankararaman2009 ; marchetti2013 . In particular Sankararaman et. al sankararaman2009 constructed dynamical equations for the concentration and the polarization field, and for the height of the film thickness. Activity was found to have two main effects on the evolution of the height field: (i) a splay induced flow that tilts the free surface and (ii) an active contribution to the effective tension. The latter can be understood by noting that active stresses pull (contractile) or push (extensile) the fluid in the direction of the long axis of the particles, giving additional elastic contribution to the stretching along that axis. By stability analysis arguments they found that, for contractile stresses, splay destabilizes the surface, while the activity contribution to tension tends to stabilize it. For extensile activity the opposite happens.

The previous results are illustrated in Figs. 9 and 10. Here the onset of instability and spontaneous flow are shown for the polar active binary mixture described in Section 2.4, in which an active polar gel (red) coexists with a passive component. Homeotropic anchoring is set both at the lower bound of the channel and at the interface between the two fluid components. The system is numerically studied by means of a hybrid Lattice Boltzmann method which combines a LB treatment for the Navier-Stokes equation with a finite difference algorithm to solve the order parameters dynamics (Section 3.5). In order to study stability with respect to bending, the interface is initially modulated by a sinusoidal perturbation. This determines a weak splay pattern in the polarization field, as shown in Fig. 9(a) and in its inset, in which white arrows represent the polarization field. As expected, the extensile system is stable under the initial splay deformations, so that, by increasing activity, these are replaced by stationary bending patterns (see Fig. 9(b)). They are accompanied by macroscopic flows, as it can be seen by looking at the velocity field, denoted with black vectors in the zooms of panels 9(b) and 9(c), with the magnitude of velocity growing linearly with $\zeta$. Then, further increasing activity, the polarization field becomes unstable (Fig. 9(d)) and the flow looses laminar character. Bends in the pattern give rise to non uniform fluxes, generating complex structures in the velocity field. A different behavior results with contractile activity. When it is strong enough, it tightens splay deformations of the initial condition of the polar field until, as shown in Fig. 10, the initial sinusoidal shape of the interface between the two fluids is completely lost and replaced by an irregular profile driven by splay polarization deformations. Fig. 10 also shows the presence of two defects of opposite charge ($+1$ and $-1$), framed with two black squares. They strongly influence the velocity pattern of the system, as shown in the inset where a quadrupolar flow can be observed in their neighborhood.

The first papers where spontaneous flow was systematically analyzed numerically are marenduzzo2007 ; marenduzzo2007pre . A slab of active nematic liquid crystal confined between two fixed parallel plates at a distance $L$ was studied by a hybrid version of LBM, with different anchoring conditions. The active nematic model is the same of that described in Section 2. The dynamics of the order parameter $\underline{\underline{\vec{Q}}}$ is governed by Eq. (18), with $\mathbf{\Xi}$ replaced by $\underline{\underline{\vec{Q}}}$ and $\mathbf{S}$ given in the last row of Table 2, with an extra active term, besides the active stress term in the Navier-Stokes equations, of the form $\lambda\underline{\underline{\vec{Q}}}$ on the right-hand side of Eq. (18). This term was suggested on the basis of symmetry considerations and also obtained by a microscopic derivation in ramaswamy2003 . Though a linear term in the nematic stress tensor also appears in the molecular field, the extra term here introduced can be regarded as active, also because no counterpart is included in the stress tensor. Positive (negative) values of $\lambda$ enhance (attenuate) self-aligning features of the nematic network, so that $\lambda>0$ can be chosen to model actomyosin suspensions at high concentration, and $\lambda<0$ to model dilute emulsions.

The main results concerning the occurrence of spontaneous flow are summarized in Fig. 11 for two different system widths $L=100\textrm{,}200$, which confirm the presence of a transition between a passive and an active phase as predicted analytically. Flow properties in the active phase are reported not to depend on the value of $\lambda$. For small $\zeta$ there is no flow and the polarization field is homogeneous. If $\zeta$ is strong enough, the system sets in the active phase, where a spontaneous flow is observed, while decreasing $L$ leads to a reduction of the active region in the parameter space. Alongside the activity parameters $\zeta$ and $\lambda$, the other key parameter is the flow alignment parameter $\xi$ (see Table 2). In fact, the transition is attained for sufficiently extensile suspensions, in the case of flow-aligning ($|\xi|>1$) liquid crystals, and for sufficiently contractile ones for flow-tumbling materials ($|\xi|<1$). In the flow-aligning case the velocity profile is characterized by the presence of bands, i.e. areas of constant shear rate, separated by narrow regions where the shear gradient reverses, similar to shear bands in non active materials fielding2004 with the number of wavelengths in the channel increasing with $\zeta$. Flow tumbling materials rearrange themselves so that only the two boundary layers flow in steady state. Simulations with periodic boundary conditions show additional instabilities, with the spontaneous flow appearing as patterns made up of convection rolls. Boundary conditions for the model in marenduzzo2007 are described in detail in marenduzzojnnm , while the numerical method in marenduzzo2011 . The phase diagram was studied, for a quasi-$1d$ system, in edwards2009 , extending previous works to the whole $(\xi,\zeta)$ plane, varying also the initial orientation of the director field.

A detailed numerical study of the dynamical spontaneous flow transition in polar active films (not by LBM but directly integrating dynamic equations) is presented in giomi2008 . In this work the effects of varying concentration were explicitly taken into account. The free-energy of the model is similar in spirit to that of Eq. (2.4) but only one phase for the concentration of the active fluid is considered (the free-energy is at most quadratic in the concentration field and no phase separation can occur). The transition to spontaneous flow is characterized by a phase diagram in a plane of two variables, related to activity $\zeta$ and to a parameter controlling self-advection777In giomi2008 , together with the active stress tensor $\sigma^{\textrm{active}}_{\alpha\beta}$ of Eq. (7), the term $\tilde{\beta}(\partial_{\beta}P_{\alpha}+\partial_{\alpha}P_{\beta})$ is also considered. This, primary to polar systems, arises from “self-propulsion” of the active units, taking into account higher order contributions in gradients in the coarse-graining onlyonlyonlysimha2002 procedure that leads to Eq. (7). This extra contribution complements the modified advected term in the evolution equations of the order parameters (discussed at the end of Section 2.5) that allows for the description of self advective phenomena.. For high values of activity and self-advection parameters a phase charaterized by spontaneous periodic oscillatory banded flows is observed. The latter, accompanied by strong concentration inhomogeneities, can also arise in active nematics, although with a physically distinct origin.

We finally mention another LBM study on spontaneous flow bonelli2016 , where the effects of a phenomenological term $\tilde{\lambda}\vec{v}$ in the dynamic equation for the polarization field $\vec{P}$ (Eq. (18)) were studied. This term is meant to favor the alignment between the polarization field and the velocity typically observed in experiments. It was shown that the interplay between alignment and activity gives rise to very different behaviors under different boundary conditions and depending on the contractile or extensile character of the fluid. With periodic boundary conditions and extensile swimmers, an almost uniform alignment between $\vec{v}$ and $\vec{P}$ is attained, while this is not anymore true for contractile systems. In the former case, when $\zeta$ and $\tilde{\lambda}$ are high enough, domains of full alignment are spaced out by small stripes or, in turn, regions in which bend distortions in the polarization field correspond to an almost zero velocity field. The behavior is yet different in bounded contractile systems. These are characterized by mutual alignment of $\vec{P}$ and $\vec{v}$, except for weak distortions close to walls for small and intermediate values of $\zeta$; moreover, fluctuations are enhanced while increasing $\tilde{\lambda}$.

The occurrence of spontaneous flow may also be accompanied by the formation of topological defects. In fact the dynamics of order parameter and velocity fields are interconnected through a feedback loop (see Fig. 10). The hydrodynamical instabilities give rise to lines of distortions in the order parameter field that are unstable to the formation of defect pairs giomi2013 . In thampi2014_1 ; thampi2014_2 an extensile active nematic has been considered and the dynamics of defects characterized. Two main stages have been identified: first, ordered regions undergo hydrodynamic instability generating lines of strong bend deformation that relax by forming oppositely charged pairs of defects. Then, annihilation of defect pairs of different charge restores nematic ordered regions which may then undergo further instabilities. In passive liquid crystals the coupling between the order parameter and the flow has significant effects on the motion of defects, generating a more intense flow around positively charged defects than for negatively charged ones Thot2002 . This phenomenon is still present in active liquid crystals, as suggested by the quadrupolar flow centered around the $+1$ defects in the inset of Fig. 10. The presence of activity gives rise to an even richer phenomenology. Full defects hydrodynamics in $2d$ polar active fluids was studied by lattice Boltzmann simulations with a hybrid scheme in Elgeti2011 . In this paper it was found that extensile activity favors spirals and vortices, like the defect highlighted by a square in Fig. 9d, while contractile activity favors aster-like defects in the polarization field like the ones boxed in Fig. 10. Defect-defect interactions have also been studied. In a contractile fluid two asters repel each other reaching a steady state with a fixed distance, that increases at larger activity. In the extensile case two asters turn into two rotating spirals, leading to a final state where the rotation continues at approximately constant rotational velocity. For low activity the angular velocity increases with $\zeta$, while above a critical value an oscillatory behavior is observed, where half clockwise rotation is followed by half anticlockwise one, resembling the previous cited oscillatory and then chaotic behaviors appearing in spontaneous flow transition for high activity888The possibility of having rotating clusters of swimmers has been also considered in the context of models of active brownian particles petrelli2018 . It has been shown that aggregates of active dumbbells, at sufficiently high activity, rotate with angular velocity inversely proportional to the average radius of the cluster.  marenduzzo2007 .

Confinement of polar/nematic pattern triggers the formation of defects, thus their dynamics is found to be particularly rich in drops of active fluids. This is not surprising given the complexity of defect topology known for passive liquid crystals rassegnadifettiliquidcrystalsgocce . So far, numerical studies available are mostly concerned with two-dimensional systems, and we will give a review in the next Section.

5 Self-propelled droplets and active emulsions

Geometrical constraints and spontaneous flow are known to significantly alter the hydrodynamics of active fluids. LB simulations show, for example, that, if a sample of active gel is sandwitched between two parallel plates, the onset of a spontaneuos flow crucially depends on the anchoring of the director field at the walls marenduzzo2007 . If the active fluid is confined in a spherical geometry, such as that of a droplet, the physics is even richer. Joanny and Ramaswamy joanny2012 , for instance, have theoretically demonstrated that the active stress can generate flows driving the spreading process of a droplet, whose shape is significantly affected by the nature of topological defects.

Recent experiments have shown that the presence of an active fluid can favour droplet self-propulsion through several mechanisms, based on chemical reactions zhang2015 ; agladze1984 , spontaneous symmetry breaking and Marangoni effects sanchez2012 ; guillamat2018 ; herminghaus2014 ; fadda2017_2 . Such experiments motivated numerous theoretical studies, with the aim of developing minimal models capable of capturing features of particular relevance in biology (as active droplets can mimick the spontaneous motion of cells bray2000 ; giomi2014 ; tjhung2012 ; blow2014 ) or in the design of bio-inspired materials sanchez2012 ; herminghaus2014 . In this section we review the theoretical studies in which the physics of an active gel under spherical confinement is described in terms of the continuum Eqs. (13)-(18)-(20), numerically solved by means of LB simulations.

Giomi et al. giomi2014 have demonstrated that those equations can describe the spontaneous division and motility of an active nematic droplet surrounded by an isotropic Newtonian fluid, a situation closely resembling the functioning of living cells. Motility results from the combination of the initial elongation of the droplet (whose stability is guaranteed by a balance between viscous and pressure forces exherted by the flow on the droplet interface and the resistance due to interfacial tension) and the instability of the active fluid to splay deformations. When the active stress is sufficiently high, splay instability dominates and triggers the formation of a spontaneous flow promoting self-propulsion. Low values of the surface tension disrupt the force balance, and lead to overstretching and eventually division of the droplet.

The generic mechanism of motility of active droplets has been also investigated analytically in Whitfield2016 and numerically in tjhung2012 by using the polar theory. Here it is shown that the self-propulsion of a contractile droplet stems from the spontaneous symmetry breaking of polarity inversion symmetry, which, in turn, triggers the formation of intense splay distortions that act as a source of kinetic energy leading to motion. The symmetry breaking can be viewed as a continuous nonequilibrium phase-transition from a non-motile to a motile state observed for a sufficiently high activity (Fig. 12). This model offers a simplified example of a cell, as a droplet containing an actomyosin solution, and suggests that motility can arise solely because of myosin contractility, rather than from its combination with actin polymerization, as often occurs bray2000 . The same hydrodynamic theory has been used to model the physics of an active polar droplet confined on a solid substrate, a situation resembling that of a crawling cell nature2015 . Here the droplet is made of an active polar fluid with contractility throughout, but actin polymerization confined in a layer close to the substrate. Such minimal description has been proven capable of capturing shapes (such as the lamellipodium bray2000 ; tjhung2015 ) andan exhibiting self-motile properties. More specifically, planar anchoring induces rotational motion in contractile droplets, while normal anchoring has the same effect in extensile ones. In both cases rotation stems from an active torque, due to a pair of bulk elastic distortions whose forma motility regimes (such as oscillatory modulations of shape bray2000 ; tjhung2015 ) by only considering a few key ingredients, i.e. actin polymerization, myosin contractility and interface anchoring.

More recently Fialho et al. fialho2017 investigated the effects of the interface anchoring of the polar field on $2d$ active polar droplets, by means of LB simulations. They found that, for large enough activity, droplets subject to strong anchoring start to spontaneously rotate, rather than exhibiting self-motile properties. More specifically, planar anchoring induces rotational motion in contractile droplets, while normal anchoring has the same effect in extensile ones. In both cases rotation stems from an active torque, due to a pair of bulk elastic distortions whose formation is controlled by a careful balance between activity and interface anchoring conditions.

As already mentioned in Section 2.3, the active gel theory has been further extended to include chirality, as this is supposed to play a relevant role in many biological systems Zhou2014 ; Naganathan2014 . Actin filaments, for example, are generally twisted in a right-ended direction Depue1965 so that, while pulling, myosin motors will tend to rotate them, creating a torque dipole. The effects of microscopic chirality on the motility of $3d$ contractile active polar droplets have been considered in Tjhung2017 . Here it has been shown that the combined effect of an (achiral) active force dipole with a (chiral) torque dipole gives rise to a rich dynamical behavior, including helical swimming, run-and-tumble motion and oscillatory swimming, the latter also observed in nature2015 . In general, while linear swimming occurs when achiral contractility is dominant, helical trajectories are observed when chiral activity becomes relevant. Such dynamic features may resemble those observed in single-celled protozoa, such as Toxoplasma gondiileung2014 .

So far we have analyzed the behavior of a single active droplet and the mechanisms leading to its self-motility properties. Active emulsions are another challenging class of systems with many potential novel applications. Examples might range from a two-fluid emulsion encapsulating active matter used for drug delivery to a biomimetic material, such as a soft tissue made up of highly-packed active droplets capable to resist to intense deformations. Active emulsions may even play a fundamental role in overcoming current major challenges affecting the design of active devices, such as controlling fuel arrival and waste removal sanchez2012 ; herminghaus2014 ; this is indeed an essential feature to create sustainable active materials, namely systems keeping their active properties over long periods of time. An active emulsion has been recently realized in the experiments by Sanchez et al. sanchez2012 , where a bundle of active nematic liquid crystal network (more specifically microtubule bundles activated through kinesin motor proteins) is squeezed at the interface of an aqueous droplet emulsified in an oil background. This paves the way towards microscopic confinement of active matter.

LB simulations would be the ideal tool to numerically investigate the hydrodynamics of a system of many active droplets and active emulsions. One may use, for instance, a multiphase approach in which each droplet is described by a different field. A different possibility is to consider negative values of $k$ in the model presented in Eq. (2.4). This models the presence of a surfactant in a binary mixture, as discussed in Section 2.4, thus favouring emulsification of the two phases. By confining polarization in one of the phases, the active behavior is then restricted to a limited portion of the system. Following this approach, the effects on morphology and hydrodynamics have been studied in a symmetric emulsion, made of an equal amount of active and passive fluid, by Bonelli et al. bone2017 , while Negro et al. negro2018 focused on highly asymmetric active emulsions, in which the active phase represents only the $10\%$ of the overall composition. Passive lamellar and hexatic orders, respectively for the symmetric and asymmetric cases, are enhanced by small contractile activity, while at more intense active dopings, an emulsion of passive droplets in an active background appears for the symmetric case. In the extensile case, a morphological transition from a stationary state to a state with chaotic patterns is mediated by the appearence of rotating droplets at intermediate active doping, regardless of the amount of active material with respect to passive phase.

Spontaneous motility in emulsions where passive droplets are immersed in an active background has also been studied in softmatter2014 , via LB simulations. The crucial mechanism leading to generation of active flows is driven by anchoring of polarization on the droplet surface. When this is homeotropic, the formation of a nearby hedgehog defect (due to the conflict anchoring with the orientation of the polar field outside the droplet) is able to drive an active flow which propels the passive droplet forward. If the gel is extensile, motility can even occur with tangential anchoring, favoured by bending deformations in the surrounding active fluid.

An important field, still in its infancy, is the study of chiral systems. Despite the motility effect of active torque has already been analyzed, as we previously discussed, intrinsically chiral active systems have gain interest only recently metselaar2018 ; carenza2019 . This route deserves a much deeper investigation: firstly because most biological extracts (such as acto-myosin systems, microtubule bundles, DNA, etc.) exhibit indeed an intrinisc cholesteric arrangement, often related to motility properties in living systems, then because the striking variety of competing metastable topologies usually observed at small pitches sec2012 ; fadda2017 is potentially useful for the design of energy-saving active soft composites built from isolated motile cholesteric droplets.

Active droplets also have the potential to exhibit a rich rheological response under shear (from a glassy to a superfluid behavior foffanoPRL ; PhysRevE.83.041910 ; cates2008 ; PhysRevLett.97.268101 ), making them suitable candidates for designing materials with new mechanical properties. A crucial question is the following: How the elastic and plastic response of an active emulsion is affected by active stress and external forcing? Decisive dynamic parameters are viscosity, shear rate, elasticity and droplet concentration. The rheological behavior of active fluids will be discussed in the next Section.

6 Rheology

The flux generated by the presence of active constituents in a suspension can deeply interact with external flows modifying the rheological response. Several experimental studies have shown, in some regime, a decrease in viscosity for pusher (Bacillus subtilis and Escherichia coli) and an increase for puller (Chlamydomonas reinhardtii) suspensions. However, a comprehensive overview of these studies is found to be far richer and more complex saintillanreview2018 .

Sokolov et al. sokolov2009 confined Bacillus subtilis bacteria in a quasi-$2d$ liquid film. Their experiments gave viscosity estimates below that of the solvent at low concentrations. At higher densities the viscosity was found to increase and exceed that of the solvent. Precise measurements using Escherichia coli bacteria were obtained by Gachelin et al. gachelin2013 and Lopez et al. PhysRevLett.115.028301 . In the first case bacteria are indeed found to decrease the viscosity below that of the solvent for low shear rates. However, the relative viscosity, given by the ratio of viscosity of the solvent in presence and in absence of suspended bacteria, increases with the shear rate, reaching a maximum above unity before shear thinning again. Lopez et al. performed their experiments in a conventional Taylor-Couette rheometer specifically build to handle low torques and viscosity. The relative viscosity displayed a trend with a low shear-rate plateau with relative viscosity less than one, followed by shear thickening. Increasing the density of bacteria was found to decrease the value of the plateau towards zero, meaning an almost vanishing value for the relative viscosity. The measured viscosity, within experimental uncertainty, actually vanishes in oxygenated suspensions for volume fractions greater than $1\%$. This surprising superfluid-like behavior should not be seen as a violation of thermodynamic principles, as the bacteria consume chemical energy. It instead suggests that the viscous dissipation in the flowing suspension is macroscopically balanced by the input of energy coming from swimming, thus allowing for a sustained flow without any applied torque. Lopez et al. also analyzed the transient stress response upon start-up and cessation of shear flow. Turning on shear, they found a sharp increase in viscosity, followed by an exponential decrease towards steady values. When shear is switched off, the shear stress undergoes a drastic reduction leading to negative values for the apparent viscosity before relaxation to a null-stress steady state.

The case of puller swimmers has been addressed by Rafai et al. PhysRevLett.104.098102 , who measured the viscosity of a suspension of Chlamydomonas reinhardtii in a Taylor-Couette flow. The suspension viscosity was always found to exceed that of the medium and to increase together with concentration. The effect of activity was spotted comparing the results obtained with suspension of dead and living swimmers; in the latter case the suspension is always found to be significantly more viscous than the one with the same volume fraction of dead cells. Weak shear thinning was also reported for high shear rates. A direct comparison of three different types of swimmers was recently performed by McDonnell et al. C4SM02742F , considering Dunaliella tertiolecta, Escherichia coli, and mouse spermatozoa. In each case they compared alive and dead cells. The viscosity measured were always above those of suspending medium, with living cells algae suspensions more viscous than dead ones, while both Escherichia coli and spermatozoa appeared less viscous when alive than dead.

The basic mechanism for viscosity modification in a suspension of microswimmers was first explained by Hatwalne et al. hatwalne2004 . It is sketched in Fig. 14. Under an applied shear profile (left panel) flow-aligning extensile systems (middle panel) enforce the applied flow, decreasing the viscosity of the suspension, while the flow generated by contractile systems opposes to the external flow thus resulting in an increase of the apparent viscosity.

They derived the coarse-grained equations governing the rheological behavior with several predictions later validated by simulations and experiments. Applying spatially uniform oscillatory shear with frequency $\omega$ in the $x-y$ plane a linear analysis hatwalne2004 gives for the $x-y$ component of the stress tensor the expression

[TABLE]

with $v_{xy}=\frac{1}{2}(\partial_{x}v_{y}+\partial_{y}v_{x})$, $L_{1}$ the nematic elastic parameter in the single constant approximation (see Table 1), $\tau_{Q}\sim 1/L_{1}$ the relaxation time of the order parameter, $\eta$ the viscosity of the suspending medium, and $f$ the active force strength. We see that viscosity, which is obtained taking the limit $\omega\rightarrow 0$, depends on the sign of $f$ and also on the geometry of the microscopic constituents, represented by the parameter $\xi$ ($\xi>0$, $\xi<0$ and $\xi=0$ for rod-like, disk-like and spherical particles, respectively). Viscosity is lowered in a suspension of rod-shaped particles by extensile activity ($f>0$), and is increased by contractile activity ($f<0$). For disc-like particles the opposite holds. Liverpool et al. PhysRevLett.97.268101 did a similar study for polar active gels, obtaining constitutive equations for the stress tensor in the isotropic phase and in phases with liquid crystalline order. The stress relaxation behavior in the various phases was discussed. Polar suspensions under uniform shear flows were further studied by Giomi et al. in PhysRevE.81.051908 , analyzing the rheological behavior in relation to the shape of the swimmers and their activity. The analytical treatment of the stress-strain linear regime demonstrates that activity lowers the linear viscosity of both extensile, rod-shaped particle and contractile, disk-shaped particle suspensions, while it increases the viscosity of contractile, rod-shaped particle and extensile, disc-like particle suspensions. For large activity the rheological response is not expected to be well described by linear analysis. Thus, Giomi et al. also integrated equations by means of Euler methods assuming uniformity in the flow direction. Depending on the value of activity and strain rate, they suggested the possibility of more exotic scenarios including hysteresis, yield-stress and non-monotonic stress-strain behaviors, and a superfluid phase with vanishing viscosity. The first two regimes have not been confirmed so far by other experiments or simulations.

In Guo201722505 a minimal phenomenological model to explain the occurrence of vanishing viscosity was proposed. This is motivated by the observation of unusual symmetric shear banding profiles in bacterial suspensions, and takes into account detailed stress balance between local viscous shear stress and active stress, and ergodic sampling of asymmetric profiles. Shear gradients in conventional complex fluids, such as worm-like micelle solutions and colloidal suspensions olmsted2008 , are always positive, while in active fluids, due to the local energy injection, one can also have regions of the system with velocity profiles characterized by negative slopes. According to the proposed model, a sheared active fluid samples all allowed shear-banding configurations, leading to a symmetric yet nonlinear shear profile.

The possibility to have states with negative viscosity, and the formation of macroscale unidirected flows is of gargantuan importance especially for applications, aiming at the design of active fluid powered motors vizsnyiczai2017 . These features have been recently explored in two papers. Slomka et al. Slomkashear2017 numerically investigated the flow pattern formation and viscosity reduction mechanisms by considering the effective model presented in Section 2 (Eq. (8)). They found, in addition to almost-vanishing viscosity phases, evidence of spontaneous formation of persistent unidirectional flow. In another paper, Loisy et al. liverpool2018 predicted negative viscosity regimes studying a model where the polarization dynamics is coupled to the flow field obeying the Stokes equation. Hydrodynamic equations are solved by finite difference methods and parameters of the model are matched with experimental data, allowing a quantitative agreement with experimental behavior.

Lattice Boltzmann methods have been also largely used to analyze numerically the behavior of active suspensions under shear. Actually, the first use of the term superfluidic for describing the behavior of an active fluid is found in the paper by Cates et al. PhysRevLett.101.068102 based on hybrid LBM simulations. Here the rheology of a slab of active gels is studied close to the isotropic-nematic (I/N) transition. For contractile gels and free boundary conditions (case in which the director is free to rotate at the boundary planes) a divergence of the apparent viscosity was found at the I/N spinodal. Close enough to the spinodal, extensile gels enter a zero-viscosity phase of N/N shear bands, with two regions of well defined nematic order and distinct shear gradients marenduzzo2007 . In the nonlinear regime, both extensile and contractile materials show non-monotonic effective flow curves, with details strongly dependent on initial and boundary conditions. Passive nematic liquid crystals ($\zeta=0$) exhibit I and N bands of high and low viscosity fielding2004 . Activity affects this behavior, widening (extensile) or narrowing (contractile) the range of shear rates for which bands are stable.

Other LBM studies on rheology of active fluids have been done by the group in Edinburgh. In PhysRevE.83.041910 , a $2d$ extensile nematic gel was studied in more detail, comparing results between planar anchoring (fixed boundary conditions) and free boundary conditions. In absence of applied shear, convection roll patterns would appear due to spontaneous flow, as discussed in Section 4. These structures, in the case of free boundary conditions and low activity, are destroyed even at low shear rates. The flow pattern becomes unsteady and spatially nonuniform. In this case shear stress takes both positive and negative values, and is characterized by large fluctuations, thus making difficult a quantitative evaluation. Increasing activity at low shear rates, the system displayed viscosity values close to zero, with less stress fluctuations. At high shear rates the laminar flow is reestablished. For very high activity the chaotic patterns observed without applied shear are only slightly perturbed at small strain rates, with the appearance of a noticeable upward curvature in the curve of the stress as a function of the shear rate. In the case of fixed boundary conditions, at low activity and shear rates the flow curves display a linear regime, while for high shear rates or high activity the results are in line with previous case. An attempt to link micro- and macro-rheology of active nematics was done in the work by Foffano et al. foffano2012 . To do so, they compared macroscopic shear numerical experiments with the results of simulations where a spherical probe particle was dragged through the active fluid. As already discussed, in sheared nematics the effective viscosity measured by bulk rheology depends on anchoring conditions and also on the system size $L$. The effective viscosity increases with $L$ in contractile fluids and decreases with $L$ in extensile ones. Normal anchoring enhances the apparent viscosity with respect to that measured at zero activity, both for contractile and extensile systems. Micro-rheology results agree with the macro-rheological characterization, and give additional informations on the role of local fluid structures on viscosity measurements. As a first important result, when a spherical particle is dragged through the active fluid, the drag force does not depend linearly on the probe radius, violating the Stokes law. When the director field is anchored tangentially to the particle surface contractility increases the drag, while extensile activity reduces it in line with the macro-rheology results. In the case of normal anchoring at the probe surface, pulling orthogonally to the far-field director leads to results similar to the tangential anchoring, while dragging the particle along the director leads the particle to move faster in contractile than in extensile nematics.

Rheology of active suspensions has been found so far to be rich of striking properties with promising perspectives towards future innovative applications. However, despite the basic mechanisms behind negative viscosity states and superfluidic regime are understood, this subject still deserves further analysis. Before such systems can be actually implemented for the design of technological applications, it is thus necessary to clarify some fundamental points such as stability of superfluidic and/or negative viscosity regimes as well as which conditions are necessary to keep the system under control. Numerical simulations, and especially lattice Boltzmann models, may play a fundamental role in further developments to this extent.

7 Active turbulence

Spontaneous flow in active suspensions evolves, for sufficient strength of the active component, into complex patterns for the velocity field that, at qualitative level, looks chaotic, resembling that of simple fluids in the passage from laminar to turbulent behavior, as suggested by streamlines of the velocity field in Fig. 3. This observation, as mentioned in the Introduction, first came from experiments conducted on highly concentrated bacterial suspensions, where the velocity field was found to develop coherent vortical structures as well as transient jets at high velocities kessler1997 ; mendelson1999 ; dombrowski2004 . Jet-like fluid motion was found in dombrowski2004 , with the speed of the jets $>100\mu m/s$, significantly larger than the speed of an individual bacterium ($\sim 10\mu m/s$). The fluid pattern observed on length scales of $100-200\mu m$ is reminiscent of a von Karman vortex street, which arises in simple fluids when the Reynolds number is much larger ($\sim 50$) compared with that of a single bacterial cell ($\sim 10^{-4})$. This justifies the terms bacterial turbulence and active turbulence, used to denote this kind of flow at low Reynolds numbers dombrowski2004 . Similar behaviors have also been reported in experiments on suspensions of both puller and pusher swimmers dombrowski2004 ; sokolov2007 ; lauga2012 , as well as microtubules bundles sanchez2012 or acto-myosin systems and microalgae drescher2010 .

A theoretical description of this complex behavior was first presented by a phenomenological model based on the Stokes equations, adapted to take into account the swimming mechanism of bacteria, each acting as a dipole stress on the fluid wolgemuth2008 . This model, solving the equations in two dimensions by using realistic parameters, empirically reproduces the observed velocity field.

The first attempt of a quantitative analysis of turbulent-like behaviors was done by Wensink et al. in wensink2012 , with a combination of experiments, simulations of self-propelled rods (SPR) and continuum theory aimed at identifying the statistical properties of self-sustained meso-scale turbulence in dense suspensions of Bacillus subtilis. Experimental and numerical data coming from SPR simulations for the energy spectrum, exhibit power-law scaling regimes better observed in two dimensions for both small ($k^{5/3}$) and large ($k^{-8/3}$) wave numbers (see Fig. 15); however, the power-law in the inertial energy range differs from the characteristic $k^{-5/3}$ decay of $2d$ high Reynolds number turbulence, as it will be discussed later. Similar results were also found making use of the continuum theory of Toner and Tu toner1998 , supplemented with the Swift-Hohenberg stress tensor presented in Eq. (8), containing higher order derivatives of the velocity gradient. Experimental data in $3d$ look qualitatively similar, but show an intermediate plateau region absent in two-dimensional systems, to indicate the spreading of kinetic energy over a wider range of scales. The same model was used in a subsequent paper dunkel2013 to consistently reproduce velocity statistics and correlations in a highly concentrated $3d$ suspension of Bacillus subtilis, while in dunkel2013_2 the linear stability analysis was presented.

Beside bacteria, eukaryotic cells with self-motile properties also show self-sustained flows. In creppy2015 a suspension of spermatozoa was found to develop a directed energy cascade characterized by a power-law scaling $k^{-3}$ at high wave numbers, a behavior that is also found in $2d$ turbulent flows and is due to non linear transfer of enstrophy, rather than energy Boffetta2012 . Recently, active turbulence features were also found in a system of self-assembled ferromagnetic Nickel microparticles dispersed at water-air interface, while subjected to an external oscillating magnetic field kokot2017 . Self-assembled spinners locally inject energy in the solvent via generation of local vortex flows, thus leading to the subsequent energy cascade towards larger scales. The hydrodynamic state is characterized by a power-law decay $k^{-5/3}$ of the energy spectrum at low packing fractions, but when this increases, steric and magnetic interactions become important so that even the exponent starts to deviate, while the system undergoes a transition toward a new phase where spinners are replaced by non-rotating clusters. The experimental observations were qualitatively reinforced through numerical simulations of disks suspended in $2d$ geometries whose dynamics was solved using the multiparticle collision approach gompper2009 .

The model of Wensink et al. wensink2012 , previously discussed, has the advantage of simplicity with respect to the active gel theory of Section 2.5. In that model the scale of variation for the velocity field is set phenomenologically, by varying the ratio between coefficients in Eq. (8), thus the explicit dynamics of the active components is neglected. Beside being pioneristic in the characterization of bacterial turbulence, whithin Wensink’s approach it is not possible to derive any information on the relations between the velocity, the active component fluctuations, and the order parameter patterns with its inherent defect dynamics. A more complete analysis of these features can be performed by considering the nematic gel theory of Section 2.5 as first done in giomi2012 . Here, decay spectra for kinetic energy and enstrophy were calculated for a relatively small $2d$ system with periodic boundary conditions. Making use of finite difference methods to solve the Beris-Edwards and the Navier-Stokes equations of Section 2.5, Giomi carefully analyzed the statistics of vortices in giomi2015 , finding that their areas are exponentially distributed and giving insights into the relations between the topological properties of the nematic director and the geometry of the flow. In particular he found that the chaotic regime is due to a feedback mechanism between the advection of nematic $\pm 1/2$ disclinations and vortical flows, whose structure is in turn closely related to presence of disclinations in their neighborhood, as demonstrated in Refs. Thot2002 ; giomi2014_2 . The number of both $+1/2$ disclinations and vortices is found to be linearly proportional to the strength of the active injection, while creation and annihilation rates of oppositely charged defects exhibit a quadratic dependence. Defects are advected along the edge of vortices at an angular velocity that is linearly dependent on the activity. As a consequence, the annihilation lifetime of disclinations is inversely proportional to the angular velocity hence, to the activity. The energy and enstrophy spectra were respectively observed to decay as $k^{-4}$ and $k^{-2}$, in the limit of high wave numbers, with no significant dependence neither on the value nor on the kind of activity (extensile or contractile).

The importance of relation between defects and vorticity has also been highlighted in Ref. putzig2016 , where a phenomenological continuum theory for overdamped active nematic liquid crystals is introduced. Three different non-equilibrium steady states are found according to the mutual effects of active torques ($\sim\nabla\times\nabla\underline{\underline{\mathbf{Q}}}$) and active convection: when the first is dominant over the second, an array of nematic defects is developed, where the comet tails of $+1/2$ disclinations are aligned throughout the system; when, instead, active convection dominates over torques, an undulated nematic free-of-defect state is formed. Nevertheless, if these contributions are equally important, a turbulent nematic state emerges, being characterized by a sharp increase both in vorticity and in the number of defects.

The growing interest in active turbulence over the last decade is motivated by some peculiar aspects: firstly, fluid active systems mostly evolve in the low Reynolds number regime, so that the development of a turbulent state due to hydrodynamic energy transfer is highly surprising; secondly, the search for an universality class for this phenomenum is extremely challenging. Moreover, even the physical mechanism that drives active systems toward the turbulent state is remarkable: energy injection takes place on small scales $l$ of the same order of the length of microswimmers, setting an upper bound for the spectral range of energy injection $k_{l}\sim 2\pi/l$, while the typical wavenumbers $k$ over which the turbulent flows take place, namely those regions where energy spectra exhibit power-law decay, are such that $k/k_{l}<1$. Hence, energy, that has been injected by the active components on microscales, is transported on much larger scales, through complex mechanisms. This is somehow reminescent of two-dimensional classical turbulence in which energy pumped into the system by an external force at a certain scale $l_{f}$ is not dissipated by viscosity, but transferred to larger scales giving rise to an inverse energy cascade Boffetta2012 . By requiring a scale independent energy-flux, one finds the Kolmogorov solution for the energy spectra $\sim k^{-5/3}$, in the inertial range $k\ll k_{f}=2\pi/l_{f}$. Yet another inertial region is expected in the direct-cascade range $k\gg k_{f}$, where the requirement of constant enstrophy flux gives a solution for the energy spectrum $\sim k^{-3}$.

In spite of the efforts taken up to now, a complete characterization for the variety of behaviors observed in the scope of bacterial turbulence is still lacking. Indeed, beside the qualitative similarity between bacterial and classical turbulent flows, it is evident, from the wide phenomenology presented before, that arguments valid for the classical high-$Re$ turbulence do not hold for active matter systems, because of the plurality of behaviors observed, nor any alternative theory has already been able to fully provide a consistent explanation of such behavior.

To fully characterize turbulent flows it is fundamental to clarify phenomena at the base of exchange of energy between different scales. For instance, as said above, in Kolmogorov theory of turbulence, energy is transported by means of the advective term leading to direct (inverse) energy cascade in $3d$ ($2d$) systems. An insight into the mechanism due to energy injection and dissipation in active fluids has been outlined by Bratanov et al. in Bratanov2015 . In this work the spectral properties of the same model as in wensink2012 are analyzed. A power-law energy spectrum was found at large scales (commonly addressed as energy-containing range), even in absence of an inertial range, namely a region of constant energy flux. The presence of further non-linearities introduced in the model, with respect to the only advective term in the Navier-Stokes Eq. (13), provides additional freedom to the energy exchange between different scales, due to forcing and energy consuming terms. Despite the model manages to replicate some phenomenological features of active turbulent systems, advective non-linear transfer plays a significantly role, that one would expect to be substantially negligible in low-Reynolds number environments. Such non-trivial hydrodynamic behavior has as consequence that exponents of the energy power-law scaling do not exhibit any universal character, since they depend on the rate and on the mean length-scale of energy injection.

Recently an insight into bacterial turbulence has been given by Doostmohammadi et al. in Doostmohammadi2017 and by Shendruk et al. Shendruk2017 , where continuum equations for the active nematic theory of Section 2.2 were solved through a hybrid lattice Boltzmann method. Shendruk found that, by confining active nematics in a microchannel, the system exhibits various morphological behaviors depending on the intensity of active doping. At small activities they found an unidirectional spontaneous flow regime, followed by an oscillating laminar flowing state, as active doping is increased, that becomes turbulent for enough intense values of the activity parameter. One of the results presented in this work is the appearance of an intermediate state between the laminar regime and the turbulent one, named Ceilidh dance. This represents the motion of paired disclinations in the nematic pattern moving along the channel, that, advected by the vortical flow, exchange partners, producing a dynamical ordered state that is reminescent of Ceilidh dancing. The transition to the turbulent state is fully analyzed in the work of Doostmohammadi et al. where vorticity puffs are used to characterize the transition. Unlike the inertial puffs that drive high-$Re$ systems toward a full turbulent state and are initiated by external forcing, in this case puffs are driven by the internal injection of energy on small scales. Far from the turbulent transition, the turbulent fraction, namely the area fraction occupied by the active puffs, is almost null, since they often vanish before splitting, but when the active doping becomes more intense while approaching the critical threshold, their lifetime grows and splitting becomes more likely. As shown in Fig. 16, if the critical threshold is exceeded the turbulent fraction grows with power-law $\sim(A-A_{cr})^{0.275\pm 0.043}$, where the adimensional active number $A=\sqrt{\zeta h^{2}/K}$, being $h$ the channel width and $A_{cr}$ the turbulent threshold. The exponent characterizing the transition closely matches the universal critical exponent of the directed percolation ($\beta=0.276$). Moreover, this is in turn very close to the exponent measured in Couette flows for inertial turbulence Lemoult2016 (since in that case $\beta=0.28\pm 0.03$). A similar characterization of $3d$ turbulence in a confined active nematics has been performed in shendruk2018 , where the crossover from quasi-$2d$ to $3d$ turbulence is controlled and driven by the deformation and twisting of disclination lines.

We finally observe that many numerical simulations performed so far Bratanov2015 ; wensink2012 ; linkmann2018 show energy transfer between different lengthscales due to hydrodynamic non-linearities, but they are not able to correctly reproduce the low-Reynolds number regime observed in experiments. Moreover active turbulence is mediated by energy injection due to self-organization of active constituents, so that one could expect that active terms and elastic interactions may play some role in energy transfer, as it happens in elastic turbulence Larson2000 ; Burghelea2006 ; morozov2007 , where the turbulent state results from elastic instabilities in viscoelastic (passive) media. No studies have been able to clarify this point yet. The works of Shendruk and Doostmohammadi and Yeomans had shed light on the route to follow to fully charachterize the turbulent-like behavior, but still a number of questions has to be answered. Is active turbulence actually turbulence? Which mechanism drives the system towards the chaotic state? How energy is exchanged between different lengthscales? Is it possible to confine and tune the turbulent-like behavior? These are matters for further research.

8 Particles and fluids

In this last section we briefly review simulation studies investigating the modeling of particle suspensions (such as colloids) dispersed in a fluid, in the context of active matter and performed via Lattice Boltzmann methods.

Inspired by natural examples such as bacteria and sperm cells marchetti2013 , a growing interest has been dedicated on producing synthetic particles capable of a similar autonomous motion ebbens2016 ; patterson2016 ; elgeti2015 ; bechinger2016 . Systems studied in experiments range from bi-metallic nanorods, self-propelled via electrophoresis (by catalytically decomposing hydrogen peroxide paxton2006 ; wang2006 ), to spherical active Janus particles, acquiring motion through self-diffusiophoresis, electrokinesis ebbens2014 and bubble nucleation wang2014 ; gregory2015 .

The numerical studies performed so far can be divided in two broad classes. The first one, on which large part of the recent research has been addressed, encompasses those made up of active particles dispersed in a passive fluid, while the second one, whose physics remains still largely unexplored, includes those in which passive particles are dispersed in an active medium.

Several models, pertaining to the first group, have been proposed to simulate microscopic swimmers interacting with a fluid at low Reynolds number. Ramachandran et al. ramachandran2006 , for instance, described a swimmer by extending the method introduced by Ladd ladd1994_1 ; ladd1994_2 ; ladd2001 , in which, for the first time, a solid particle was computationally modeled within a LB mesh. In this approach a solid particle can be created by initially defining a closed surface, represented by a set of boundary links joining neighbouring fluid nodes, located outisde and inside the surface. The former are fluid nodes, whereas the latter are solid nodes. At the surface, stick boundary conditions, bouncing back the density of the fluid moving along a boundary link, are implemented. For computational convenience the interior of the particle remains fluid (i.e. lattice nodes inside and outside the particle are treated in the same way), a condition overall acceptable for a single-phase fluid as inertial effects are negligible999Extensions of the model, in which the internal fluid is excluded, have been also implemented nguyen2002 . . The active contribution is included by adding a force, exerted by the particle to the fluid, distributed on each boundary node located between the fluid and the solid particle. The sum of these forces must be zero for the Newton’s third law, a condition fulfilled solely by a force dipole. Self-propulsion is then achieved by modifying the symmetry of these forces. An asymmetric force distribution generates a net force on the particle (due to the difference of the fluid velocity in its back and in its front) and promotes motion, whereas a symmetric distribution does not yield a force unbalance, although the particle still generates a non-zero quadropolar velocity field in the fluid. The former particles are called movers and the latters are called shakers.

A complementary mechanism promoting particle self-propulsion was proposed in Ref. Llopis2006 , and consisted in adding a constant amount of momentum with fixed magnitude to the particle. To restore momentum conservation, the same amount of momentum is subtracted to the fluid nodes connected to the solid ones. This approach has been used to study the collective dynamics of self-propelled particles dispersed in a fluid solvent, and was found to reproduce, for instance, the formation of transient aggregates of particles, as well as the transition of the particle mean square displacement from ballistic towards diffusive motion at low concentration Llopis2006 . A similar mechanism has been also used to investigate the hydrodynamics of active rotors Llopis2008 , systems capturing, for example, the rotating motion of molecular motors on the cell membrane noji1997 .

An alternative approach has been afterwards used to simulate suspensions of swimming particles (again built using the Ladd’s method) starting from the squirmer model introduced by Lighthill lighthill1952 ; blake1971 , in which particle’s motion is now triggered by a predefined axis-symmetric tangential velocity distribution, imposed on the surface of the particle alarcon2013 . By using such velocity distribution as boundary condition of the Stokes and the continuity equations (as the inertia of the fluid is neglected), the mean fluid flow induced by a minimal squirmer can be computed as

[TABLE]

where the terms on the right hand side derive from an expansion with Legendre polynomials truncated to the second mode. Here $\theta$ is a polar angle and $\theta=0$ defines the swimming direction of the squirmers. The dimensionless parameter $\epsilon=B_{2}/B_{1}$ defines the type of squirmer: $\epsilon>0$ describes a puller (or contractile squirmer) while $\epsilon<0$ describes a pusher (or extensile squirmer). The case with $\epsilon=0$ corresponds to an apolar squirmer (or shaker).

Lattice Boltzmann has been also applied to investigate the hydrodynamics of active three-bead linear swimmers earl2007 , a system for which an analytic solution is available najafi2004 . Here the swimmer is constituted of three spheres of predefined radius linked with sufficiently thin rods to neglect hydrodynamic interactions. Lengths and angles between the rods can be modified in a periodic and time irreversible manner by the action of internal forces and torques, which favor change in the swimmer shape and potentially motion when coupled to the fluid. In this model there is no net interface separating the fluid from the bead, whose interior, as in the model proposed by Ramachandran et al. ramachandran2006 , is then essentially fluid. The swimmer-fluid interaction is incorporated in three steps, the first of which computes the total linear and angular momenta of the swimmer. Afterwards position and shape of the swimmer are updated with respect to the original one (such that linear and angular momentum are conserved) and finally the motion of the swimmer is coupled to the fluid to obtain updated fluid velocities used to calculate the equilibrium distribution functions of the LBM earl2007 ; pooley2008 . Such approach was found to reproduce analytical results with sufficient accuracy, in particular when swimmers with a high number of beads (even more than three) and close to each other (or near a surface) are considered. An alternative and perhaps computationally simpler way to model solid object on a grid was suggested by Smith and Denniston smith2007 . Here the surface of the swimmer is represented by a large number of point particles, which ensure a non-slip condition by introducing a drag force between the particles and the nodes of the lattice Boltzmann mesh pooley2008 . Although this method was found efficient in simulating the swimming velocity of a single swimmer, it suffered of large discretization errors affecting inter-swimmer hydrodynamic interactions.

Due to a limited computational efficiency, the models described so far have been rarely applied to simulate highly anisotropic particles. This drawback has been recently circumvented by de Graaf et al. degraaf2016 ; degraaf2016_2 , that modeled active colloids of arbitrary shape as clusters of spheres coupled to the LB fluid through the scheme introduced by Ahlrichs and Düenweg ahlrichs1999 . Unlike the grid-based Ladd’s method, in this approach a particle is described through a set of points coupled to the fluid through a frictional force, which depends on the relative velocity (in analogy to the Stokes formula for a sphere in a viscous fluid) between the fluid and the points. Due to this coupling and to the interpolation of the force on the LB mesh at the particle position, a hydrodynamic shell forms around the points, which now acquire an effective hydrodynamic radius. A solid object is then obtained when a high enough number of points is considered degraaf2015 ; degraaf2015_2 . Self-propulsion is achieved by means of a persistence force applied along a direction vector assigned to the particle, while an equal and opposite counter force is applied to the fluid nash2008 ; hernandez2005 ; saintillan2007 . Its location sets the nature of the swimmer (whether contractile or extensile) and the structure of the fluid flow in its surrounding. Besides being a facile approach, this method has the advantage to incorporate hydrodynamic interactions with higher order multipole moments, (in addition to the usual dipole ones), and has been found to well reproduce far-field theoretical results in system with periodic boundary conditions and in a spherical cavity with no-slip walls deGraaf2017 . A likewise computationally efficient method, in which swimmers are described using the point-force implementation developed by Nash et al. nash2008 , has been used to perform unprecedented large-scale LB simulations of $\sim 10^{6}$ hydrodynamically interacting swimmers in a cubic box with periodic boundary conditions stenhammar2017 . In this work the authors show that swimmers move in a correlated fashion well below the transition to bacterial turbulence, and elaborated a novel kinetic approach capturing such behavior with results in quantitative agreement with LB simulations.

More recently, lattice Boltzmann methods have been extended to simulate more complicated physical systems in which active particles play, once more, a relevant role. A remarkable example is the process of cross-streamline migration (often called margination schmid1980 ) of stiff active particles (such as synthetic nanoparticles used for drug-delivery) migrating towards the vessel walls when moving in blood flows. Such phenomenon, studied in Ref. gekle2016 , is attributed to hydrodynamic interactions of red blood cells with stiff particles, and disappears when red blood cells are absent. Active particles are modeled by triangulated spheres whose internal grid nodes are massless points flowing with the local fluid velocity and connected each other through stiff harmonic springs and bending potential to preserve the grid arrangement. This approach is often referred as immersed boundary method vahid2014 ; vahid2015 . A force is then applied to the center of mass of the particles (which become active) and along the opposite direction with respect to the fluid flow.

A much less investigated system is that in which a passive particle is embedded in active fluid. This has been done in Refs. foffanoPRL ; foffano2012 , where LB simulations show a violation of the Stokes’ law when a particle (modeled by using the Ladd’s method) is dragged in an active fluid. Even more strikingly, a negative viscous drag in the steady state of a contractile fluid is found for large enough particles. Such simulated microrheological experiment, in principle realizable in the laboratory, highlights the fact that, although at its infancy, the study of these systems may unveil novel and intriguing physics, potentially useful for the design of novel active soft materials.

Clearly many efforts have been addressed to model systems which (i) correctly describe the hydrodynamics of active particles and (ii) capture new dynamical and mechanical properties of great importance for future practical applications. Several directions of research can be envisaged starting from the results achieved so far. A largely unexplored field is that of the active rheology, in which suspensions of particles (e.g. passive) are subject to an external perturbation (such as a shear flow) in an active liquid crystal. Here the nature of the active system (either contractile or extensile) as well as the particle volume fraction and the shear stress may dramatically affect the rheological response. In addition, the presence of topological defects, due to the conflicting anchoring of the liquid crystal with that on particle surface, may foster the formation of turbulent-like velocity patterns, whose physics is still under investigation. LB simulations performed so far usually neglect thermal fluctuations, which may however be crucial especially when the particle size decreases up to nanometers. Even more intriguing seems the possibility of a multiscale coupling of the LB with other simultation methods, such as dissipative particle dynamics or immersed boundary method, to capture, for instance, the near-contact interactions of active colloids at the interface of a binary fluid.

9 Conclusions

In this review we have offered an introduction to theoretical approaches in continuum modeling of active fluids together with the description of the most appropriate lattice Boltzmann techniques for numerical simulations of this new fascinating category of soft matter. This has been done in Sections 2 and 3, respectively. One of the most immediate consequences of activity is the occurrence of spontaneous flow. We have illustrated in Section 4 that this emerges as a result of bending instability of the polarization field in extensile systems and of splay deformations in contractile fluids. The role played by nematic or polarization defects has been also discussed. The presence of geometrical constraints and/or external forcing makes the behavior of active fluids or mixtures even richer, with many interesting and somehow unexpected features. Active fluids confined in droplets can spontaneously flow or spread faster, if compared with similar passive cases, and this opens the way to many futuristic applications coming from the design of new bio-inspired materials. Self-propelled droplets can be even model systems for living cells and this also justifies the enormous recent research interest towards this matter, reviewed in Section 5. In Section 6 the rheological behavior of active fluids is discussed. Unidirected motion and superfluidic behavior are examples of results in this field whose theoretical and applicative consequences are matter for future research. We have discussed these issues with the more recent perspectives of research in these fields. Large self-propulsive forces induce large flows with chaotic velocity patterns. This behavior has been termed active turbulence even if the possibility of a quantitative description similar to that of standard turbulence at high Reynolds number is still a matter to be properly understood. The status of research in this field has been reviewed in Section 7. Finally, in Section 8, it is shown how mixed LBM-particle algorithms can be used to simulate the effects of the active bath on particles in suspension. In all these subjects numerical simulations have demonstrated to be an essential tool for the characterization and for the general understanding of the phenomena considered. Between different numerical methods the role of LBM has been prominent as shown by the considerable amount of results in literature coming from LBM studies. The reason of this is mainly in the easy adaptation of lattice Boltzmann algorithms to the study of complex fluids and active fluids. We hope that this review can be useful to other researchers in this new field of soft matter.

Acknowledgments

We warmly thank L. Biferale, P. Digregorio, D. Marenduzzo, E. Orlandini, I. Petrelli for many fruitful discussions. Simulations supporting the theories presented in this review were run at Bari ReCaS (https://www.recas-bari.it/) e-Infrastructure funded by MIUR through PON Research and Competitiveness 2007-2013 Call 254 Action I, at ARCHER UK National Supercomputing Service (http://www.archer.ac.uk) through the program HPC-Europa3 and at Marconi (CINECA - http://www.hpc.cineca.it/) under CINECA-INFN agreement (Project No. INF19-fldturb).

Contribution Statement

All the authors have equally contributed to conceive and write this review. L.N.C. and G.N. have performed the simulations needed for Figs. 3,6,7,9,10.

Appendix A Biaxial nematics

In this Appendix we briefly discuss how to extend the uniaxial order parameter to biaxial nematics and how biaxiality provides information about the localization of topological defects in nematic liquid crystals.

A biaxial nematic is a nematic liquid crystal with three distinct optical axis and, unlike an uniaxial liquid crystal, it does not have any axis of complete rotational symmetry. Hence one can define three perpendicular axes ${\bf n}$, ${\bf m}$ and ${\bf l}$ (two are sufficient, since the third one would be perpendicular to the others), or director fields, for which there is a reflection symmetry. The order parameter in three dimensions is then

[TABLE]

where $S$ and $\varepsilon$ are called scalar order parameters. This representation of $\underline{\underline{\vec{Q}}}$ is a traceless symmetric second order rank tensor with five independent components. If the smaller of the two scalar order parameters is very small (like in many systems), one recovers the $\underline{\underline{\vec{Q}}}$ tensor in the uniaxial approximation.

For a biaxial nematic the three eigenvalues are in general different, and the diagonal representation of the $\underline{\underline{\vec{Q}}}$ tensor is

[TABLE]

where $\varepsilon\leqslant S$, with $\varepsilon$ gauging the degree of biaxiality.

The investigation of nematic defects in the context of a Landau-de Gennes theory has shown that their core presents a heavy degree of biaxiality [272]. By following the approach of Ref. [273], this can be measured by computing three scalars, $c_{l}=\tilde{\lambda}_{1}-\tilde{\lambda}_{2}$, $c_{p}=2(\tilde{\lambda}_{2}-\tilde{\lambda}_{3})$ and $c_{s}=3\tilde{\lambda}_{3}$, where $\tilde{\lambda}_{1}$, $\tilde{\lambda}_{2}$ and $\tilde{\lambda}_{3}$ (with $\tilde{\lambda}_{1}\geq\tilde{\lambda}_{2}\geq\tilde{\lambda}_{3}$) are three eigenvalues of the diagonalised matrix $G_{\alpha\beta}=Q_{\alpha\beta}+\delta_{\alpha\beta}/3$. These parameters fulfill the following properties: $0\leq c_{l},c_{p},c_{s}\leq 1$ and $c_{l}+c_{p}+c_{s}=1$. An ordered unixial nematic will give $c_{l}\simeq 1$, while the isotropic state, where both $S$ and $\eta$ are approximatively null, corresponds to $c_{s}\simeq 1$. Finally the biaxial case implies $c_{p}\simeq 1$.

Appendix B Chapman-Enskog expansion

In Section 3.2 we presented a LB algorithm to solve the hydrodynamics of a simple fluid. We show here that Eqs.(20) and (13) can be recovered in the continuum limit, starting from the evolution equation (21) for the distribution functions $f_{i}$. Two approaches can be followed. The first one starts from a Taylor expansion of the left-hand side of Eq. (21) [274], whereas the second one, discussed below, uses a Chapman-Enskog method, that is an expansion of the distribution functions about equilibrium, which assumes that successive derivatives are of increasingly high order in the Knudsen number $\epsilon=\lambda/L$. This is a dimensionless number defined as the ratio between the molecular mean free path $\lambda$ and a characteristic length $L$ of the system. This number determines whether a macroscopic continuum mechanics or a microscopic statistical mechanics formulation of fluid dynamics should be used. For small values of $\epsilon$ ($\epsilon\ll 1$) the mean free path is much smaller than $L$ and a continuum theory is a good approximation. To take into account both ballistic and diffusive scales, spatial density fluctuations of order $\mathcal{O}(\epsilon^{-1})$ are assumed to relax over time scales of order $\mathcal{O}(\epsilon^{-2})$. A suitable expansion for temporal and spatial derivatives as well as for distribution functions is

[TABLE]

built assuming that there is a diffusion time scale $t_{2}$ slower than the convection one $t_{1}$.

We start by expanding the left-hand side of equation (21) to the second order in $\Delta t$:

[TABLE]

where Eq. (23) has been used to express the collision operator. By substituting Eq. (70), (71) and (72) into Eq. (73) one obtains

[TABLE]

By retaining at most terms of second order in $\epsilon$, the previous equation reads

[TABLE]

Finally, grouping terms of same order in $\epsilon$, we get

[TABLE]

In the following paragraphs we will use these relations to recover continuum equations up to second order in the Knudsen number.

Recover Continuity Equation

To recover the continuity equation one can start by summing Eq. (76) over lattice velocities with the constraints given in Eqs. (26) and (27). One then gets

[TABLE]

and, by using again Eq. (27)

[TABLE]

By performing a summation over lattice velocities in Eq. (77), one gets

[TABLE]

which is the continuity equation at first order in the Knudsen number. To recover the complete time derivative according to Eq. (71), we need to explicitly compute the term $\partial_{t_{2}}\rho$. By applying the differential operators $\epsilon\partial_{t_{1}}$ and $\epsilon\xi_{i\beta}\partial_{\beta_{1}}$ to Eq. (77) we obtain

[TABLE]

[TABLE]

and, summing both equations:

[TABLE]

Note that the left-hand side of this equation is exactly the term in round brackets of Eq. (78), that now becomes

[TABLE]

By summing over lattice directions and using Eqs. (79) and (80), we get

[TABLE]

that, together with Eq. (81), reads

[TABLE]

Finally, after restoring the canonical differential operators (through Eqs. (71) and (72)), we get the continuity equation

[TABLE]

at second order in the Knudsen number.

Recover Navier-Stokes Equations

The procedure to recover the Navier-Stokes equation is analogous, albeit less straightforward, than that used for the continuity equation. We will proceed by calculating the first-order moment of Eq. (77) and Eq. (78). First multiply by $\xi_{i\beta}$ both members of Eq. (77) and sum over index $i$, to get

[TABLE]

To get the Navier-Stokes equation to second order in the Knudsen number we need to calculate the first-order moment of equation (78). We can then multiply Eq. (85) by $\xi_{i\gamma}$ to obtain

[TABLE]

and, by summing over lattice velocities, we are left with

[TABLE]

Now we must determine an expression for the summation in square brackets. From Eqs. (76) and (77) we note that

[TABLE]

where we have used Eq. (29) in the second equality. The second term of the second line of Eq. (92) can be written in terms of the equilibrium distribution functions given in Eq. (34) and of the related coefficients in Eq. (31)

[TABLE]

while the first round bracket in the second line of Eq. (92) can be written by means of Eq. (81) and Eq. (89) as

[TABLE]

In the last line terms of order $v^{3}$ were neglected, an approximation valid as far as the Mach number is kept small. Now substituting Eqs. (94) and (93) into Eq. (92) we find, after some algebra, that

[TABLE]

This term, in turn, enters Eq. (91), which now reads

[TABLE]

Finally, summing this equation with Eq. (89) and using the canonical differential operators (i.e. Eqs. (71) and (72)), we obtain the Navier-Stokes equation

[TABLE]

where $p=(c^{2}/3)\rho$ is the isotropic pressure and the shear viscosity is given by

[TABLE]

Appendix C Recovering continuum equations for the algorithm described in Section 3.3.2

In this Appendix we show the calculations to recover the continuum equations obtained by means of the forcing method algorithm discussed in Section 3.3.2.

We start with a Chapman-Enskog expansion of the distribution functions and of the derivatives according to relations (70)-(72), and for the forcing term we assume that

[TABLE]

. By Taylor expanding the evolution Eq. (21), we get

[TABLE]

We now substitute Eqs. (70)-(72) and (100) in Eq. (101) and, after grouping terms of the same order in $\epsilon$, we get

[TABLE]

From these equations one gets the zeroth-order moments of the distrubution functions

[TABLE]

the first-order ones

[TABLE]

and the zeroth, first and second moments of the forcing term

[TABLE]

These relations can be explicitly calculated by using Eq. (45). Now combining the zeroth-order moment of Eq. (103) with Eqs. (105) and (109), one gets

[TABLE]

the continuity equation at first order in the Knudsen number. To recover it at second order we apply the differential operators $\partial_{t_{1}}$ and $\xi_{i\alpha}\partial_{\alpha_{1}}$ to Eq. (103) and then, by perfoming the difference between the two resulting equations, one has

[TABLE]

Now we substitute the latter into Eq. (104) and, by using Eq. (103), we obtain

[TABLE]

Finally summing this one over index $i$, by means of Eqs. (102)-(110), we get

[TABLE]

By summing this equation with Eq. (112) we obtain the continuity equation at second order in the Knudsen number.

To reproduce the Navier-Stokes equation, we start by computing the first moment of Eq. (103) that, after using Eqs. (105)-(111), reads

[TABLE]

Following the same procedure for Eq. (114) we get

[TABLE]

Finally, summing these two equations, one can restore the Navier-Stokes equation

[TABLE]

where we require that

[TABLE]

and the kinematic viscosity is

[TABLE]

Appendix D Boundary conditions

In many practical situations, such as in a system under shear flow, one may be interested in studying the physics of the system within a confined geometry. Here we describe the implementation of boundary conditions of a sheared bidimensional fluid defined on a lattice of size $L_{x}\times L_{y}$ and confined between two parallel flat walls located at $y=0$ and $y=L_{y}$. Two key requirements are necessary for a correct description of the physics:

• •

no flux accross the walls,

• •

fixed velocity $v_{x}^{*}$ along the walls,

which correspond to the following relations on the wall sites:

[TABLE]

Assuming a $d2Q9$ lattice geometry (see Fig. 5) with the walls located along the lattice links (i.e. along the lattice vectors $\vec{\xi}_{1}$,$\vec{\xi}_{3}$), one can explicitly write the previous relations at $y=0$ (the bottom wall):

[TABLE]

Note that after the propagation step, functions $f_{0}$, $f_{1}$, $f_{4}$, $f_{7}$ and $f_{8}$ are known, so that one can use relations (122) and (123) to determine the three unknown distribution functions $f_{2},f_{5},f_{6}$. This system of equations can be closed by adding the bounce-back rule:

[TABLE]

The two remaining distribution functions $f_{5}$ and $f_{6}$ are then given by [275]

[TABLE]

With this choice for inward-pointing distributions, the desired momentum at the boundary is achieved. Unfortunately this scheme does not allow for the local conservation of mass since, after the collision step, inward-pointing distributions are not streamed. In [276] an improvement of this scheme was proposed to overcome such a problem. In the following we will use notation $f_{pi}$ to identify the outgoing distribution function in a wall lattice site at time $t-\Delta t$, while $f_{i}$ denotes those streamed from neighboring sites at time $t$. Besides conditions in (121) it is required that the fraction of mass moving towards the wall or eventualy still on a wall site at time $t-\Delta t$ is the same that moves from the wall or stay still on the walls at time $t$. This is expressed for a bottom-wall site by the following relation:

[TABLE]

where $f_{0}$ must be determined by solving the system of Eqs. (121) and (127) together with the bounce-back condition (124). This leaves unchanged the solutions for the unknown $f_{5}$ and $f_{6}$ in Eqs. (125) and (126), but provides a new expression for $f_{0}$ that is thus given by:

[TABLE]

Such scheme can be easily adjusted to the case of the pure forcing method presented in Section 3.3.2. The only difference lies in the momentum conservation relations [277] that in such case read as follows,

[TABLE]

The system of Eqs. (129) together with Eq. (127) admits the following solutions:

[TABLE]

where the outward-pointing distribution $f_{2}$ was fixed by the bounce back condition (121).