Mesoscopic simulations at the physicschemistrybiology interface
Abstract
We discuss the Lattice BoltzmannParticle Dynamics (LBPD) multiscale paradigm for the simulation of complex states of flowing matter at the interface between Physics, Chemistry and Biology.
In particular, we describe current largescale LBPD simulations of biopolymer translocation across cellular membranes, molecular transport in ion channels and amyloid aggregation in cells. We also provide prospects for future LBPD explorations in the direction of cellular organization, the direct simulation of full biological organelles, all the way to up to physiological scales of potential relevance to future precisionmedicine applications, such as the accurate description of homeostatic processes.
It is argued that, with the advent of Exascale computing, the mesoscale physics approach advocated in this paper, may come to age in the next decade and open up new exciting perspectives for physicsbased computational medicine.
Contents
 I INTRODUCTION
 II BOLTZMANN KINETIC THEORY
 III LATTICE BOLTZMANN FOR CONTINUUM HYDRODYNAMICS
 IV LATTICE BOLTZMANN FOR GENERALIZED HYDRODYNAMICS
 V FLUIDS AND PARTICLES: THE LATTICE BOLTZMANN  PARTICLE DYNAMICS SCHEME
 VI SIMULATIONS AT THE PHYSICSCHEMISTRYBIOLOGY INTERFACE
 VII SIMULATIONS AT THE PHYSICS/CHEMISTRY/BIOLOGY INTERFACE: DREAMING AHEAD
 VIII PCB Modeling versus Big Data Science
 IX SUMMARY AND PERSPECTIVE
 X Acknowledgements
 X.1 High Performance Simulations of Particle Dynamics
 X.2 Achieving High Performance for Lattice Boltzmann methods
 X.3 Overlap between computation and communication
 X.4 Sparse and irregular geometries
I Introduction
Thanks to the spectacular advances of computer technology (hardware and software) on the one side, and mathematical modelling on the other, in the last few decades modern science has come to the point of providing a quantitative description of many biological systems, whose complexity would have been regarded as missionimpossible until only recently. In this Review, we shall illustrate the point through several concrete examples.
Notwithstanding such major advances, the challenge of modelling biological and physiological systems remains formidable, as it actually amounts to cover some ten decades in space (from molecules to the human body) and easily twice as many in time. No mathematical/computational model in the foreseeable future can take up such a challenge headon, i.e. by direct simulation of all the actual mechanisms, scales and levels involved in the process.
Coarsegrained methods come in many flavours and families, depending on the range of scales and problems they are targeted to, but in this review we shall focus on a specific mesoscale technique, known as the Lattice Boltzmann method (LBM), namely a minimal lattice version of the Boltzmann equation Boltzmann (2012); Cercignani and Berman (1976) which has witnessed a burgeoning growth for the description of complex flow phenomena across an impressively broad range of scales.
The Lattice Boltzmann Equation (LBE) was developed as a computational alternative to the discretization of the NavierStokes (NS) equations of continuum fluid mechanics Benzi et al. (1992); Succi (2001).
Over the years, however, it has made proof of an amazing and largely unanticipated versatility and ability to describe a broad variety of phenomena involving complex states of flowing matter, beyond the strict realm of continuum hydrodynamics, including nontrivial flows at micro and nanoscale. Thanks to this versatility, and to the coupling with various families of mesoscale particle methods, in the last decade LBE has gained increased status for the simulation of many complex flow problems at the interface between fluid dynamics, chemistry, material science and biology. These include, for instance, multiphase and multicomponent flows with complex interfaces, the motion of suspended bodies under strong geometrical confinement, possibly with chemical reactions Succi (2018); Krüger et al. (2017). After revisiting the main ideas behind the Lattice Boltzmann (LB) theory, in this review we discuss current and future prospects of multiscale/level Lattice BoltzmannParticle Dynamics (LBPD) simulations at the physicschemistrybiology interface, in an attempt to identify and portray outstanding problems of potential relevance to clinical applications in a notsodistant future, i.e. computational medicine.
Computerassisted medicine is a consolidated branch of modern science, which generally develops around two separate pillars; molecular biology and macroscale physiology. The former is heavily leaning on bioinformatics and data science tools, while the latter usually relies upon the methods of continuum and fluid mechanics. In this Review, we portray a third, alternative approach, based on mesoscale physics, i.e. physicsinformed coarsegrained models of microscale biological processes, possibly embedded within their physiological environment. This mesoscale approach is grounded into the intermediate level of the description of matter, namely kinetic theory, in both its versions, Boltzmann’s kinetic theory and Langevin stochastic particle dynamics. As a result, the main tools of the mesoscale approach are fluids, particles and probability distribution functions. The multiscale LBPD concept is thus apparent; Boltzmann naturally connects upwards to continuum hydrodynamics, while Langevin connects downwards to the molecular level. Combining the two in a single computational harness, opens up a direct route from the continuum to the atomistic word, and back.
Several specific examples have already shown the potential of LBPD simulations in areas straddling across Physics and Biology. A selected set of applications will be described in detail in this Review, to provide a taste for the breadth of applications that have been tackled in the recent past. Readers keen on specific details are kindly directed to the original literature.
Clearly, the PCB interface is enormously rich and varied, ranging from nanometric macromolecular phenomena, to peptidic aggregation and biopolymer translocation, to cellular motion and active matter, as sketched in Fig. 1. In a few words, the topic goes way beyond the scope of the present work. Nevertheless, it is important for the reader to appreciate the great flexibility of the method to cope with multiple scales of motion and, more importantly, its ability to incorporate the desired degree of biochemical specificity inherent to the problem at hand.
The LBPD framework can be enriched in the direction of describing complex flowing states of matter at micro and nanoscales. Although not “rigorous”, such variants prove capable of providing new physical insights into highly complex states of flowing matter. Remarkably, such extensions can be put in place without compromising the outstanding parallel amenability to parallel processing of the method. On the assumption that computing power will keep growing in the next few decades, if only perhaps at subMoore’s paces, this state of affairs spawns tremendous opportunities to gain new insights into a series of fundamental problems dealing with complex states of flowing matter in general, and with special focus on those relevant to biology and medicine.
The paper is organised in three main parts. In the first one, we discuss the basic aspects of Boltzmann’s kinetic theory with special emphasis on its lattice version for fluid dynamics and its extensions to softmatter and biological applications, including the coupling to particle dynamics for the motion of suspended bodies.
In the second part, we describe selected applications to biological and physiological systems, such as biopolymer translocation, ion channels, protein diffusion and amyloid aggregation in cellular environments. For a quick visual appreciation, see figures 2, 3, 4, to be commented in detail later in this manuscript.
In the third part, we provide a prospective view of a series of problems at the physicschemistrybiology interface, which may become accessible once Exascale computers are with us. Special attention is paid on their potential import for clinical applications, such as the direct simulation of biological organelles and the quantitative description of hemostatic processes.
Finally, due to the crucial role played by highperformance computing in this story, in the Appendix we provide an extended account of the main issues involved with the implementation of the LBPD scheme on highend parallel computers in the Exascale range.
The main message we wish to convey in this Review is that a mesoscale physicsbased approach to computational medicine may come to age in the next decade.
Ii Boltzmann Kinetic Theory
The Boltzmann equation (BE) is the core of Boltzmann’s kinetic theory, that, in turn, is the cornerstone of nonequilibrium statistical mechanics, a pillar of theoretical physics at large Boltzmann (2012). Besides its paramount conceptual value as a bridge between the microscopic world of atoms and molecules and the macroscopic world of thermohydrodynamic fields, the BE also provides a concrete tool for the quantitative investigation of a broad variety of practical nonequilibrium transport problems Cercignani and Berman (1976). However, the BE is all but an easy piece to work with: a nonlinear integrodifferential equation in (phasespace plus time) dimensions. This motivates a relentless search for new methods to solve the BE either analytically or numerically, the latter option usually covering a broader ground. Graeme Bird’s Direct Simulation Monte Carlo (DSMC) method has played a leading role in this respect and continues to do so to the present days Bird (1994). In principle, DSMC solves the BE directly and in full, i.e. accounting for the specificity of molecular interactions, as well as strong nonequilibrium effects, using a stochastic particle technique, whence the Monte Carlo label. This comes at a major computational cost, which is why various approximations have been developed and considerably refined over the years Ketsdever and Struchtrup (2016); Dimarco et al. (2018). Close to local equilibrium and away from confining elements, however, molecular details become increasingly irrelevant. Universality takes stage and more economical descriptions can be devised. The basic idea is to relinquish the “irrelevant” details while still preserving the basic properties of macroscopic physics, namely the symmetries and conservation laws which secure the emergence of the NS equations from the underlying molecular dynamics. Among others, a description which has gained major interest for the last three decades is the LB method Benzi et al. (1992); Succi (2001, 2018); Krüger et al. (2017). LB was devised with the specific intent of providing an alternative to the discretization of the NS equations for the numerical solution of continuum hydrodynamic problems. This still is its mainstay and, for some authors, also the only place where it belongs. Indeed, the use of LB for flows beyond NS was ruled out Luo (2004); Junk et al. (2005), mostly on account of the lack of a rigorous asymptotic limit. The above nogo has been proven largely overrestrictive and nowadays applications beyond the strict realm of continuum fluid dynamics abound, especially in the direction of soft matter. Since problems in biology and medicine hardly involve fluid mechanics alone, such developments are of direct relevance to computational explorations at the interface between physics, chemistry and biology, the main scope of this Review.
ii.1 The Boltzmann Equation
The Boltzmann equation (BE) of classic kinetic theory is basically a continuity equation in sixdimensional phasespace, namely Boltzmann (2012):
(1) 
where is the probability density of finding a molecule at position in ordinary space, with velocity at time . The lefthand side of the BE represents the streaming of the molecules under the effect of a force field and reflects the Newtonian mechanics and of a representative molecule, while the righthand side describes intermolecular interactions. For the case of a dilute gas, as originally considered by Boltzmann, these interactions are typically taken in the form of twobody local collisions, since higher order encounters are much less frequent, hence negligible in the collision count. Even under such major simplifications, solving the BE presents a very daunting challenge on account of its highdimensionality, six phasespace dimensions plus time, as well as due to the nonlinear (quadratic) integral character of the collision operator Cercignani and Berman (1976).
Regardless of the complexity of the underlying microscopic interactions, the collision operator must comply with the massmomentumenergy conservation laws, namely:
(2) 
In addition, it must also secure compliance with the Second Principle, which amounts to supporting a socalled Htheorem, namely:
(3) 
In other words, the dynamics of the distribution function must converge to a universal global attractor, corresponding to the thermodynamic equilibrium.
The macroscopic fluid variables are obtained by a linear and local contraction of the Boltzmann distribution, namely:
(4)  
(5)  
(6) 
where is the mass density, the flow speed, and the flow temperature in spatial dimensions.
Central to the emergence of hydrodynamic behavior is the notion of local equilibrium. This is defined as the specific form attained by the Boltzmann distribution once collisions come in complete balance, i.e.
(7) 
Inspection of the Boltzmann collision operator provides the following universal MaxwellBoltzmann (MB) local equilibrium distribution:
(8) 
where is a normalization constant in spatial dimensions, is the thermal speed and
is the peculiar speed, i.e., the molecular velocity relative to the fluid one, in units of the thermal speed. The reader familiar with statistical mechanics will readily recognize the canonical distribution in the comoving frame of the fluid, with the identification .
A few comments are in order.
First, the MB distribution depends on space and time only through the hydrodynamic fields, , its dependence on the velocity variable being a universal Gaussian distribution. This is a strict consequence of Eq. 2, i.e., the microscopic conservation laws.
Such dependence is largely arbitrary, with the caveat that it should be weak on the molecular scale. More precisely, the macrofields should not show appreciable changes on the scale of the molecular meanfree path, that is
(9) 
where designates any macrofield and is the molecular mean free path. The above ratio, known as Knudsen number, serves as the smallness parameter controlling the emergence of the hydrodynamic limit from Boltzmann’s kinetic equation. Ordinary fluids dynamics holds in the range and below.
ii.2 From Boltzmann to NavierStokes hydrodynamics
The conceptual path from BE to the NS equations of continuum fluids is based on two fundamental steps:

Projection of the Boltzmann equation upon a suitable basis function in velocity space, typically Hermite polynomials in cartesian coordinates,

Multiscale expansion using the Knudsen number as a smallness parameter, on the assumption of weak departure from local equilibrum.
The projection generates a hierarchy of partial differential equations for the kinetic moments:
(10)  
(11)  
(12) 
where
(13) 
and denotes the th order tensor Hermite polynomial. Note that is a tensor of rank , namely (scalar) is the fluid density, (vector) is the fluid current and (second order tensor) is the momentumflux, whose trace delivers (twice) the kinetic energy of the fluid and the triple tensor is the flux of momentum flux.
The left hand side clearly shows that the moment hierarchy is open, since the time derivative of is driven by the divergence of . This is the mechanism by which heterogeneities fuel nonequilibrium. Also to be noted that the righthandside of the first two equations is zero because collisions conserve mass and momentum. However, they do not conserve momentumflux, which is why the righthandside of the third equation is nonzero, expressing the relaxation of the momentumflux to its equilibrium expression. Such relaxation takes place on a collisional timescale , which in turn fixes the kinematic viscosity of the fluid. Instantaneous relaxation () denotes the infinitely strong collisional regime whereby collisions do not leave any chance to nonequilibrium to survive, formally corresponding to the idealized case of a perfect (zero dissipation) fluid.
Macroscopically, this corresponds to the inviscid Euler equations.
All real fluids, though, relax in a short but finite time (strictly speaking this is also true for superfluids), and consequently the moment equations present an open hierarchy which needs to be closed, somehow. This is where the assumption of weak departure from local equilibrium takes stage.
By that assumption, one formally expands the distribution function and spacetime derivatives in powers of the Knudsen number, replaces the expansion in the moment equations, and collects homologue terms order by order in the Knudsen number. To zero order, the Euler equations are obtained, whereas the first order delivers the NS equations of dissipative fluids, namely
(14)  
(15) 
where is a second order tensor formed by spatial derivatives of the flow field. Typically:
(16) 
where is the symmetrized gradient tensor, is the divergence of the flow and is the unit tensor. The first scalar is the dynamic shear viscosity whereas associates to the bulk viscosity, .
Despite their deceivingly simple physical content, essentially mass and momentum conservation (Newton’s law), as applied to a finite volume of fluid, the NS equations prove exceedingly difficult to solve, as they involve the nonlinear evolution of a threedimensional vector field, often in a complex geometry set up. This sets a formidable challenge to even the most advanced computational methods, whence the ceaseless hunt for more efficient numerical methods. Those methods include a broad array of techniques to discretize the NS equations, using finite differences/elements/volume schemes.
Three decades ago, however, an entirely different route was devised, which consists in attacking fluid dynamics “from the bottom”, i.e., appealing to a microscopic description of the fluid states matter, namely a highly stylized version of molecular dynamics known as Lattice Gas Cellular Automata Chopard and Droz (1998); Rivet and Boon (2005). In a nutshell, the idea is to introduce a boolean lattice fluid, consisting of a set of boolean particles whose dynamics is confined to the lattice sites. Boolean, here, means that the state of the system at a given lattice site and instant in time is uniquely defined by a set of binary digits, coding for the absence/presence of a corresponding particle moving with unit speed along one of the links connecting each lattice site to its neighbours. By a suitable choice of the lattice connectivity and interaction rules, such Boolean system can be shown to reproduce the NS equations of continuum fluid dynamics. An associated LB equation was also derived in the process of taking the boolean automaton to NS, but its computational capabilities went under noticed. Even though the LGCA did not make it into a competitive tool for computational fluid dynamics, it nonetheless paved the way to the idea of computing fluid flows by simulating fictitious particle dynamics instead of discretising the NS equations. LB Higuera et al. (1989) fully inscribes within this line of thought.
ii.3 Boltzmann equation for Biology?
The Boltzmann factor is a household name in biology, as it governs the statistical behaviour of a broad class of equilibrium and nonequilibrium (activated processes) phenomena of utmost relevance to biological systems. But, how about the Boltzmann equation? At first sight, Boltzmann kinetic theory, in its original form at least, offers little scope for biological applications, since it formally applies to dilute states of matter in which molecular collisions are rare, the socalled weakly coupled regime, in which manybody interactions can safely be neglected with respect to binary molecular encounters. On the contrary, most biological phenomena are hosted primarily by condensed and soft matter systems in which manybody effects play a primary role.
Here, however, a change of perspective proves exceedingly fruitful. Rather than the original Boltzmann equation for actual molecules, we shall refer to model Boltzmann’s equations for fluidparticles, the latter denoting the effective degrees of freedom describing the behaviour of representative groups of molecules. Historically, the distinctive feature of model Boltzmann equations is a dramatic simplification of the collision operator, in an attempt to relinquish most mathematical complexities while retaining the essential physics at hand. The most popular Boltzmann model is the celebrated BhatnagarGrossKrook equation, in which the collision operator is replaced by a simple singletime relaxation term Bhatnagar et al. (1954):
(17) 
where is the local equilibrium and is a relaxation frequency controlling the relaxation to the local equilibrium on a timescale .
The key advantage of this change is neat: model equations are more flexible by construction, hence they can be modified (extended) not only to simplify the collision operator but also to describe a wide host of physical effects not included in the original Boltzmann equation.
In particular, they can reinstate the effects of, manybody interactions, via effective onebody forces, in the spirit of Density Functional Theory (VlasovBoltzmann Equation) Hansen and McDonald (1990) statistical fluctuations, through appropriate stochastic sources (Fluctuating Boltzmann Equation) Ladd (1993) farfromequilibrium inhomogeneities, via suitably extended collisionrelaxation operators in which the relaxation time is promoted to the status of a selfconsistent dynamic field Chen et al. (2003); Higuera et al. (1989); d’Humières (1992).
The three extensions above inscribe to the general framework of Reverse Kinetic Theory (RKT), the strategy whereby the kinetic equation is designed topdown, based on prescriptions securing compliance with macroscopic hydrodynamics in the first place, and then adding “molecular” details “ondemand” by the specific application under investigation.
RKT reverses the canonical bottomup route, whereby kinetic equations are derived from the underlying microscopic models and proves quite effective in bringing Boltzmannlike equations within the realm of condensed and soft matter physics. However, care must be exercised in securing compliance of the topdown approach with the basic principles of statistical physics, namely symmetries/conservation laws as well as evolutionary constraints (the Second Principle and its local form, known as Htheorem).
This all said and done, a practical question still remains: although simplified, the model BE’s still leave in an unwieldy dimensional space. Here, a timehonoured ally of any computational scientist, the lattice, makes its glorious entry.
ii.4 Hydrodynamics for Biology
Hydrodynamics and biology couple across an amazingly broad spectrum of scales, ranging from the macromolecular level in the compartments of living cells, all the way up to industrial bioreactors, rivers, lakes and oceans. From macromolecular motion in crowded cellular spaces, to the deformations of membranes, to the motion of cells and the selfpropulsion of bacteria, the common goal is to capture the effect of hydrodynamics under both equilibrium and nonequilibrium conditions. The latter may arise under the effect of an external flow, such as the transport of proteins and cells in the blood stream, or via the consumption of energy in cellular metabolic pathways. Indeed, the response of biological matter to temperature, pH or mechanical forces, plays a key role in most biological processes. For instance, biological macromolecules, cells or tissues, are usually fragile and easily damaged by hydrodynamic or shear forces, as they occur far from equilibrium and under strong confinement. Following in the footsteps of hydrodynamic experimental techniques, customarily used to measure basic molecular properties, such as weight, size and shape, computer simulations can be used in biophysical chemistry to explore and assess the mechanical and dynamical properties of macromolecules far from equilibrium, hence much closer to the in vivo conditions relevant to medical purposes. The benefits of including hydrodynamic forces is even more apparent whenever thermodynamic forces, stemming from solutesolvent interactions, can also be taken into account. A comprehensive computational framework including both hydrodynamic and thermodynamic forces, stands therefore as a most desirable target. Owing to its inherently intermediate nature, between atomistic and continuum descriptions of matter, (lattice) kinetic theory sits at a vantage point to meet this goal.
Iii Lattice Boltzmann for Continuum Hydrodynamics
The basic idea of the LBM is to constrain the velocity degreesoffreedom to a discrete lattice with sufficient symmetry to protect the conservation laws which secure the emergence of standard hydrodynamic behavior in the macroscopic limit. LB is based on the idea of representing fluid populations on a uniform, cartesian grid.
The standard LB scheme in singlerelaxation time (BGK) form reads as follows Higuera et al. (1989); Qian et al. (1992):
(18) 
where is the discrete Boltzmann distribution associated with the discrete velocity , running over the discrete lattice, to be detailed shortly. In the above, is a relaxation parameter controlling the fluid viscosity and is the lattice local equilibrium, basically the local MaxwellBoltzmann distribution truncated to the second order in the Mach number. The truncation is not a luxury, but a requirement dictated by Galilean invariance (GI).
Let us remind that GI refers to the invariance of the NS equations under an arbitrary change of the local fluid velocity ; upon such change, the NS equations stay the same, provided is replaced by . In continuum kinetic theory, Galilean invariance is encoded within the dependence of the local MaxwellBoltzmann distribution on the relative velocity of the molecules with respect to the fluid one, namely . As a result, an observer in the comoving frame, that is a frame moving at the local fluid velocity, experiences the same local equilibrium as if there were no fluid motion. To be noted that, in the above, is an arbitrary function of space and time, indicating that Galilean invariance is a local and continuum symmetry, i.e., it holds even if different regions of the fluid move at different velocities, which is of course the case in most fluids of practical interest. Galilean invariance is reflected by the specific form taken by the moments of the equilibrium distribution, and specifically by those explicitly relevant to hydrodynamics, namely mass, momentum and the momentum fluxtensor, that is:
where latin subscripts run over spatial coordinates . The above expression follows straight from the property of gaussian integrals in velocity space.
One might naively expect that the same would be true in the lattice, provided the local MaxwellBoltzmann distribution is retained with the plain replacement . Straightforward algebra shows that the situation is different not for a mere mathematical accident, but as a consequence of the fact that a local and continuum symmetry cannot remain unbroken in a discrete lattice. More precisely, it cannot remain unbroken for any arbitrary velocity field. It turns out, though, that this is possible whenever the fluid velocity is much smaller than the sound speed, i.e., in the lowMach number limit.
Under such a perturbative approximation, replacing velocity integrals with discrete sums returns exactly the same moments as in the continuum, namely:
(19) 
provided that lattice tensors up to fourth order are isotropic. Of course, this by no means implies that GI is fully restored, but simply that the GIbreaking terms are confined to kinetic moments higher than order four.
In other words, the hydrodynamic constraints can be matched perturbatively, by expanding the local Maxwellian to second order in the Mach number, which is sufficient to recover the (isothermal) NS equations, since those equations are quadratic in the fluid velocity. Full Galilean invariance for an arbitrary flow field implies instead an infinite series in the Mach number, corresponding to the full expansion of the local Maxwellian in Hermite polynomials.
The actual expression of the discrete local equilibria reads as follows:
(20) 
where and represent the dipole and quadrupole contributions, respectively. Hereafter, repeated indices are summed upon. In the above, is a set of latticedependent weights, normalized to unity, which represent the lattice analogue of the global (noflow) MaxwellBoltzmann distribution. Finally is the lattice sound speed.
We hasten to note that, at variance with its continuum counterpart, the expression (20), being a polynomial truncation of the MaxwellBoltzmann equilibria, is nonnegative definite only in a finite range of fluid velocities, typically of the order of . This configures the LB method as an appropriate description of quasiincompressible, low Machnumber flows.
The discrete velocities matching isotropy constraints up to fourth order can be shown to be subsets of the D3Q27 (DQ is a widely used notation to indicate a LB scheme in dimensions using a set of velocities) mother lattice with discrete speeds in three spatial dimensions. D2Q9 and D3Q27 are the direct tensor product of the elementary onedimensional D1Q3 stencil, , in two and three spatial dimensions, respectively (see Fig. 6).
The expansion (20) implies that hydrodynamic LB flows are bound to be quasiincompressible, i.e., Mach number well below unity. Likewise, third order kinetic moments, describing energy and heat flux, are not correctly reproduced, since these terms require sixthorder isotropic lattice tensors. Such constraints can not be met by lattices confined to the first Brillouin region described by D3Q27: higher order lattices extending beyond the first Brillouin cell are necessary. So much for the low Mach number approximation, which is specific to the lattice. But what about the low Knudsen limit, which, on the contrary, lies at the very roots of the convergence of Boltzmann to NS?
Since the Knudsen number controls the heterogeneitydriven departures from local equilibrium, it is intuitively clear that the lowKnudsen hydrodynamic limit implies further constraints on the nonequilibrium component of the momentumflux tensor, which amounts to recovering the continuum expression of the stress tensor.
Based on (19) and (20), it can be shown that in order for such constraints to be met, isotropic tensors, again up to order four, need to be exactly reproduced in the lattice. It may come as a surprise that nonequilibrium constraints can be matched at the same order of isotropy of the equilibrium ones. The reason is parityinvariance, namely the fact that each discrete lattice velocity comes with a mirror partner . Indeed, non equilibrium constraints involve lattice corrections driven by fifthorder lattice tensors, which are identically zero due to parity invariance. Readers interested in the straightforward but laborious algebraic details, may look up the vaste literature on the subjectChopard and Droz (1998); Rivet and Boon (2005); Hénon (1987).
Next, let us comment on the source term, , at the righthand side of the LB equation (18). This term stands for a generic source of mass/momentum/energy, describing the coupling of the fluid to the surrounding environment. Mass sources are typically associated with the presence of chemical reactions, turning one species into another in multicomponent versions of the LB for reactive flows.
In the case of inert flows, the source term typically encodes the momentum exchange due to the coupling to external (or internal) fields, such as gravity or more complex interactions, like selfconsistent forces reflecting potential energy interactions within the fluid, as well as thermal fluctuations. As we shall see, the two latter cases are crucial for the extension of LB beyond NS hydrodynamics.
From the operational standpoint, the source term acts like a bias promoting the populations which move along the local force field and setting a penalty on those that move against it. It is therefore clear that the strength of such term is subject to stringent stability and positivity constraints.
iii.1 From Lattice Boltzmann to continuum hydrodynamics
The conceptual path taking LB to NS is exactly the same as in the continuum theory, with the crucial caveat of turning around the (many) catches associated with lattice discreteness. To make a long story short, it amounts to securing the proper symmetries of the lattice tensors entering the set of (lattice) moment equations. Like always with lattice physics, the name of the game is to erase the lattice dependence to the desired order. For the case of isothermal, incompressible fluids, the desired order is the fourth one. In equations, and using coordinate notation for tensors Hénon (1987):
(21)  
(22)  
(23) 
where latin indices run over spatial coordinates. In the above, is a set of weights normalized to unity and oddrank tensors are automatically equal to zero by parity invariance, i.e., every discrete velocity comes with an equal and opposite partner, so that the sum of the two gives a null. Finally is a lattice dependent constant expressing the sound speed in the lattice, as per the expressions above.
With the symmetries secured, everything proceeds like in the continuum theory, with another important caveat though, namely the fact that the effective meanfree path of the lattice fluid is replaced by the lattice spacing whenever the latter is larger than the physical one, the typical case in most macroscopic LB simulations.
This is a very technical and somehow thorny issue, whose details are outofthescope of the present Review. Nevertheless, we wish to caution the reader that the physics taking place at the scale of a few lattice spacings should always be inspected with a big pinch of salt, because it is constantly in danger of breaking the hydrodynamic assumptions.
Once these catches are disposed of, one ends up with a lattice fluid obeying an ideal equation of state and a kinematic viscosity given by:
(24) 
where is the natural lattice viscosity. Note that the stability range of the discrete timemarching, , also secures the positivity of the kinematic viscosity.
Both ends of this range must be handled with care. In the lowviscosity regime , typical of turbulence, strong gradients may develop posing a serious threat to the numerical stability of the scheme. A powerful variant of the basic LB includes a selfconsistent tuning of the relaxation parameter so as to ensure compliance with local entropy growth (Htheorem) Karlin et al. (1999); Succi et al. (2002). That variant, known as Entropic Lattice Boltzmann (ELB) is normally intended to simulate highReynolds macroscopic flows, but lately is proving very effective also to stabilize microscale LB simulations with sharp interfaces.
In the opposite highviscous regime, , the one most relevant to biological applications, the viscosity may formally diverge, signalling a departure from hydrodynamic behaviour, due to the onset of ballistic motion. As a result, whenever possible, LB simulations should be kept away from both limits, say .
We shall return to this important point in the Section ”Future Challenges”.
iii.2 Boundary conditions
In the early days, boundary conditions were hailed as one of the main assets of LB, and, to a certain extent, they still are. As a matter of fact, since information always travels along straight lines, even complex geometries can be handled by comparatively straightforward computational methods based on elementary mechanical operations. For instance, noslip on solid walls can be implemented through a simple bounceback between distributions propagating along opposite directions, i.e., from the fluid to the wall and viceversa (See Fig. 7).
For the sake of generality, let us consider a fluid flowing in a bounded domain confined by a surrounding boundary . The problem of formulating boundary conditions within the LBE formalism consists in finding an appropriate relation expressing the incoming (unknown) populations as a function of the outgoing (known) ones . Outgoing populations at a boundary site are defined by the condition
where is the outward normal to the boundary cell centered in . Incoming populations are defined by the opposite sign of the inequality. In mathematical terms, this relationship translates into a linear integral equation
(25) 
where the kernel of the boundary operator generally extends over a finite range of values inside the fluid domain. This boundary operator reflects the interaction among the fluid molecules and the molecules in the solid wall. Consistently with this molecular picture, boundary conditions can be viewed as a special (sometimes even simplified) type of collisions between fluid and solid molecules. Physical fidelity can make the boundary kernel quite complicated, which is generally not the idea with LBE. Instead, one usually looks for expressions minimizing the mathematical burden without compromising the essential physics.
In particular, one seeks minimal kernels fulfilling the desired constraints on the macroscopic variables (density, speed, temperature and possibly the associated fluxes as well) at the boundary sites . This may lead to a mathematically underdetermined problem, more unknowns than constraints, opening up an appealing opportunity to accommodate more interface physics within the formulation of the boundary conditions. However, it also calls for some caution to guard against mathematical illposedness. By now, a vast literature which, however, goes beyond the scope of the present Review, is available on this topic Krüger et al. (2017); Zou and He (1997); Chen et al. (1996).
Here, we just wish to point out that the main issue controlling the complexity of the boundary problem is whether the collection of boundary points lies on a surface aligned with the grid (socalled staircase wall boundary) or it is given by a set of offgrid surface elements (sometimes called surfels). The latter case is significantly more complex and requires extracare, the immediate benefit being that the nearwall physics is secondorder accurate as compared to the firstorder accuracy of the staircased approximation (for a detailed account see Succi (2018); Krüger et al. (2017)). Advanced applications of the offlattice boundary method prove capable of dealing with fairly complex geometries, an important but highly technical topic which also goes beyond the scope of the present Review. For full details see the recent books Krüger et al. (2017); Succi (2018)
iii.3 The bright sides of Lattice Boltzmann
Why does LB represent an appealing scheme for simulating complex states of flowing matter ?
Several features merit highlight, but essentially they all track down to the benefits of working with extradimensions opened up by the sixdimensional phasespace inhabited by kinetic theory.
More in detail, the upshots are the following.
First, the information always travels along straight lines, regardless of the spacetime complexity of the emergent hydrodynamics.
What we mean by this is that the discrete distributions move in sync along the discrete lightcones, defined by:
no matter the complexity of their spacetime dependence.
Since the discrete velocities do not depend on space and time, the streaming is exact, no information is lost in hopping from one lattice site to another, literally an errorfree operation. This stands in sharp contrast with the selfadvection term of macroscopic hydrodynamics, whereby the velocity fields, or any other field for that matter, moves along spacetime dependent material lines defined by the fluid velocity itself, typically a complicated function of space and time.
Such independence explains the outstanding LB amenability to parallel computing, a practical asset which can not be underestimated.
Second, space and time always come on the same (firstorder) footing. In particular, this means that dissipation is not expressed by secondorder spatial derivatives (Laplace operator) but simply as a local relaxation to a local equilibrium, as described earlier. This is a substantial advantage for confined flows, whose dynamics is largely dictated by the presence of solid boundaries, where spatial gradients are usually at their largest. Taking space derivatives near the boundaries is notoriously prone to numerical inaccuracies, especially whenever the geometry of the domain is irregular, as it is often the case for biological flows. LB equipped with suitable formulations for curved boundaries can significantly mitigate such difficulties.
Third, LB carries pressure as just any other macroscopic field, with no need of solving a (usually very expensive) Poisson problem to compute the pressure field consistent with an incompressible flow. This is because the discrete distributions carry the momentumflux tensor “on their back” and consequently the pressure obeys its own dynamic equation. More specifically, the second order momentumflux tensor, , obeys the relaxation equation
(26) 
where is the thirdorder energy flux tensor, latin indices running over spatial dimensions. In the limit where the tensor is enslaved to its equilibrium value, the time derivative can be dropped and the energyflux tensor can be approximated by its equilibrium expression. Under such conditions, the above equation delivers .
Once the due lattice symmetries are fulfilled, the equilibrium component delivers the advection and pressure terms of the NS equation, while the third order nonequilibrium term delivers the dissipative term.
Fourth, coupling of the fluid to a broad variety of other physical phenomena, is readily achieved by formulating suitable expressions of the source term .
Typically this reflects the action of mesoscale forces describing the effective interactions between fluidfluid and fluidsolid molecules. The positive side effect is that an entire new world of complexity may be supported simply by inserting a few tens of lines of additional code.
This feature is crucial to the extension of LB beyond the traditional realm of dilute gases and particularly to biological flows, as we shall detail in the sequel.
Once again, this rosy picture only refers to the conceptual fresco, actual implementations requiring great care in dribbling lattice artefacts.
iii.4 The dark side of the LB moon
Following a witty line by Daan Frenkel, it is often more instructive to analyse what can go wrong in a computer simulation, rather than indulging in the glorious description of its success stories. To abide by this wise attitude, in the following we mention things that can still go wrong with LB simulations, the dark side of the LB moon.
LB draws much of its computational simplicity to the features of uniform and regular lattices, that, at times, recalls an ideal Legoland. Realistically complex geometries may sometimes challenge this setup, and call for more flexible and elaborate formulations, such as interpolation procedures for curved boundaries, local gridrefinement, mergers with finitevolumes techniques, to name just a few of some popular options. It is only fair to acknowledge that implementing such extensions maybe laborious and possibly tax the computational simplicity as well.
Another limitation is the finite compressibility. While buying major computational savings through dispensing with the Poisson problem for pressure, a finite amount of compressibility must be tolerated in return. Such effects are usually negligible for lowReynolds flows, but they need nevertheless to be watched carefully in open flows, especially at outlets, where spurious density waves may eventually backpropagate into the fluid and undermine the accuracy and sometimes even the stability of the simulation.
A similar observation applies to flows with strong thermal transport. For one, it should be appreciated that LB is essentially athermal, in the sense that the discrete Boltzmann distribution is represented by a linear superposition of monochromatic beams with no dispersion in velocity space, hence, zero temperature in a strict equilibrium thermodynamic sense. Nevertheless, a kinetic temperature can still be defined as a measure of the dispersion around the mean flow velocity, namely, in spatial dimensions, . The inclusion of strong thermal effects typically requires higher order lattices, accommodating sixth order tensors describing the flux of energy Krüger et al. (2017); Nie et al. (2008). Several alternatives are also possible which go beyond the scope of the present work (for details see chapter 22 in Succi (2018)). Despite major progress, it is fair to say that thermal LB schemes still lag behind their athermal counterparts in terms of numerical robustness. Consequently, LB is often coupled to independent thermal solvers typically based on finitedifference or finitevolume methods.
These limitations off the breast, we believe it is fair to surmise that the appealing features of the LB method largely overweigh its weaknesses.
Iv Lattice Boltzmann for Generalized Hydrodynamics
The ideas and methods presented so far deal with flows of “simple” fluids that can be described in terms of the NS equations (note that simple fluids can give rise to highly complex flows, turbulence being a prime example in point).
Modern hightech applications, nanoengineering and biology in the first place, set a pressing demand of quantitative understanding of more complex states of flowing matter, in which fluids interact with external or selfconsistent fields, undergo chemical reactions, phasetransitions, or interact with a variety of suspended bodies, such as colloids or biological molecules.
This emerging sector of modern science, often referred to as “complex fluids”, or more trendily, “soft matter physics”, portrays a multidisciplinary scenario whereby fluid dynamics makes contact with other disciplines, primarily chemistry, material sciences and biology as well.
There is a growing evidence that LB, and extensions thereof, holds a vantage point as a computational framework for the simulation of complex states of flowing matter. Ideally, LB would fill the gap between fluid dynamics and molecular dynamics, namely the huge and allimportant region where fluid dynamics breaks down and molecular dynamics is not yet ready to take over for lack of efficient algorithms and computing power (see Fig. 8) Boon and Yip (1991).
LB is a natural candidate to fill that gap because of its mesoscopic ability to incorporate microscopic details into the kinetic theory formalism via suitable external fields and/or equivalent generalizations of local hydrodynamic equilibria: a microscope for fluid mechanics, a telescope for molecular dynamics Succi (2018).
The extension of LB to generalised flows is based on a number of major upgrades of the basic LB theory. In the sequel, we focus on the following selection:

Reactive systems

Charged flows

Flows far from equilibrium

Fluctuating hydrodynamics

Nonideal fluids

Flows with suspended bodies
We now proceed to a more detailed discussion of the items above.
iv.1 AdvectionDiffusionReaction systems
Chemical reactivity is an essential element to deal with when facing the task of simulating systems at the PCB interface Boon et al. (1996); Succi (2018); Coveney et al. (2016a); Alowayyed et al. (2017).
Biochemical reactions lie at the heart of most biological phenomena, controlling species interconversion. They involve the breaking and making of covalent bonds, often catalyzed by enzymes, ubiquitous in all metabolic and synthetic reactions within the cell and in the body. Reactions occur in bulk conditions in singlephase (homogeneous) or multiphase (heterogeneous) environment, typically at the interface between phases or in a porouslike environment. Biological reactions take place under a wide host of different conditions. The antigenantibody binding in solution or the binding of small molecules to plasma proteins in blood, clotting reactions on the surface of blood vessels, the hydrolysis of adenosine triphosphate (ATP\nomenclature[ADR]ATPadenosine triphosphate) on the surface of endothelium, ion transport, and the release of nitric oxide are examples of homogeneous reactions. Heterogeneous reactions within tissues include aerobic metabolism, glucose consumption, and receptormediated endocytosis.
In all cases, reactions can occur either in noflow or under flow conditions, as for example the enzyme reactions on the surface of the blood vessels. Diffusion and convection affect the rates of homogeneous and heterogeneous reactions. Because reactants must diffuse or convect to the surface where they react, the mass transfer mechanisms occur in sequence with the reaction process. In many cases, the effects of reaction and diffusion on the reaction rates cannot be easily separated and numerical methods provide the only viable route to study the process.
LB is a powerful framework to capture diffusion and reaction since its mathematical apparatus is by no means confined to fluid equations.
As a matter of fact, by suitable tuning of the local equilibria, relaxation matrix and external forces, it permits to generate a very broad variety of linear and nonlinear partial differential equations, including those describing relativistic and nonrelativistic quantum phenomena Succi (2002a); Mendoza et al. (2010).
Of relevance are advectiondiffusionreaction (ADR) \nomenclature[ADR]ADRAdvectionDiffusionReaction equations of the general form:
(27) 
where is a species concentration, say molecules per unit mass or volume. The lefthand side describes the passive transport along the fluid material lines, whereas the righthandside describes diffusion plus the effects of chemical reactions occurring in the fluid moving at the barycentric velocity .
Although thermodynamically favorable, many reactions are limited by the energy barrier which needs to be crossed in order to form the activated state, so that, in the absence of the catalyst, reactions simply do not occur. In its presence, the rate of reaction can increase dramatically, although the change in energy between reactants and products is not affected. In fact, the enzymes affect the rate of a reaction, not its equilibrium. For many reactions involving a single substrate, the rate of consumption of the substrate follows the MichaelisMenten equation
where is the magnitude of the rate of disappearance of the substrate (or reactant) and is the Michaelis constant, the concentration of the substrate at which the reaction attains half of its maximum rate.
Another reaction term particularly relevant to populations biology is the logistic expression
where describes the (Malthusian) growth rate and characterizes the nonlinear decay of the species due to competition for, say, space and/or food. The ADR class falls naturally within the formalism, by simply defining a LB distribution for the concentration, such that and specifying the following local equilibria:
where both fluid and lattice speeds are normalised by the sound speed . To be noted that in the above is the fluid speed, which does not match the current , because ADR’s conserve mass but not momentum. The LB is then particularly well suited to solve the ADR equations in complex geometries, such as those occurring in morphogenesis Ponce Dawson et al. (1993); Ayodele et al. (2011), heterogeneous catalysis Falcucci et al. (2016) and related phenomena.
iv.2 Charged Fluids
Electrostatics plays a vital role in biological processes and requires handling electrolytic solutions and charged solutes typically in flow conditions, the realm of the socalled electrokinetics. Solutes range from strongly charged molecules, such as DNA, to weakly charged macromolecules, such as proteins, where partially screened mobile charges and charged surfaces are ubiquitous conditions of biological settings. Prominent examples include i) the stability of proteins as a function of pH and ionic strength, such as saltingin and saltingout effects due to the interplay of protein charges with the aqueous/saline environment; ii) proteinligand association processes, including enzymatic allostery modulated by saltbridges, saltbridges in virus assemblies and thermal stability; iii) membrane proteins and their specific electrostatic properties, whereby the collective charge of the intracellular residues tend to be more positive as compared to the extracellular ones (the socalled “positive inside” rule), and so on. The associated biological flows are driven by electroosmosis and electrophoresis. Wellknown examples of those processes include i) ion channels and the gating of ions across the cell membrane that regulate electrical signalling in secretory and epithelial cells, as much as the cell volume; ii) the electrokinetic flow due to the glycocalyx layer covering cells, a polyelectrolyte exhibiting a surface charge and responsible for an electrochemical gradient that regulates mechanotransduction; iii) the mitochondria and the chloraplast, where proton gradients generate a chemiosmotic potential, also known as a proton motive force, for the synthesis of ATP.
Although the general principles of electrokinetics are fairly well understood Masliyah and Bhattacharjee (2006), methods that translate these principles into accurate numerical predictions are still in their infancy. In fact, deriving the interactions between charged solutes and the solvent requires computing the interactions of a large number of molecules and the averaging of these over many solvent configurations. Such daunting computational requirements are partly overcome through the continuum mesoscopic approach, whereby solvent and ions are described in the continuum and in a preaveraged sense Capuani et al. (2004); Reboux et al. (2006); Wang and Kang (2010).
The LB framework solves complex electrokinetic problems through an efficient formulation of the electrolytic solution as a multispecies problem: one species for the neutral solvent, water, and two for the positively and negatively charged ionic components. The Poisson equation provides the solution for electrostatics and the selfconsistent forces for the transport of the fluid species, in the socalled VlasovPoisson approximation.
Let us briefly survey the LB method for charged multicomponent systems. The fluid mixture is composed by three sets of populations labelled by index , two ionic components with charges , being the proton charge, density and velocity . Given the barycentric velocity , the species relative velocity is . The aqueous medium accommodates solute biomolecules, with the ith particle having position and valence . The electrostatic potential is obtained by solving the Poisson equation,
(28) 
where the medium has dielectric permittivity , duly complemented by boundary conditions of the Dirichlet or Neumann kind. For insulating confining walls, where the surface has local charge density , this reads , where is the unit vector normal to the surface.
The dynamics of each species follows the evolution equation:
(30) 
and the force term
The local selfconsistent forces are , where is the drag force exerted on species , resulting in a crossdiffusion coefficient , where is a relaxation frequency.
So much for the governing equations. However, a hallmark of charged systems is that local electrostatic forces can be dramatically intense, often exhibiting rapid spatial modulations of the electrolytic densities, as in the presence of double layers.
Such occurrence can lead to severe numerical instabilities in methods based on the direct solution of NS equation, as a consequence of local violations of the CourantFriedrichLevy stability condition Courant et al. (1928). Thanks to its inherently small time step (in macroscale units), LB can often handle stiff forces without loosing stability. However, strong polyelectrolytes are problematic to handle and cannot be treated directly, thus some sort of charge rescaling is required Datar et al. (2017). Despite these liabilities, LB simulations of electrolytic systems show stable behaviour over a wide range of solute charges and molarities, in particular in the submolar range that covers a large portion of biological conditions.
iv.3 Flows far from equilibrium
Most biological systems operate far from equilibrium, i.e., they draw energy from their environment and dissipate heat back to it, thereby lowering their own entropy content at the expense of the environment. This is the operating principle of the so called “dissipative structures”, a cornerstone of nonequilibrium thermodynamics Prigogine (2017). In the process of dumping entropy to the environment, they manage to migrate across a sequence of different nonequilibrium steadystates (NESS), in which they deliver different functions.
Examples of NESS in biology include molecular machines, cells in motion, metabolic pathways and many others. For instance, the kinesin protein walks along the microtubule, carrying cargos from one part of the cell to another, by absorbing energy from ATP hydrolysis and converting chemical energy into mechanical work, of which is used for motion and the rest is dissipated to the surroundings. Offequilibrium conditions imply that biological agents exchange mass, momentum, energy and entropy with the environment through transient or steady currents and fluxes. Clearly, such complex spacetime dependent network of currents and fluxes, typically operating on a broad spectrum of concurrent space and time scales, needs to be captured by the numerical approach Takahashi et al. (2005). Whence a major need of multiscale methods.
Indeed, these currents and fluxes typically occur through thin (atomic scale) interfaces, which resolve the tension between competing mechanisms, say chemical reactions and molecular diffusion, through a sudden spatial transition between distinct bulk phases, say the liquid and vapour in a multiphase flow.
For the case of a diffusionreaction system, the width of the transition region can be estimated as
where is the diffusion coefficient and is a typical chemical timescale. By expressing , with being the meanfree path and the collisional timescale, we obtain . This shows that, unless chemical reactions are much slower than inert collisions (slowchemistry), the interface width is comparable with the molecular mean free path, or shorter (fastchemistry).
Another example in point are foams and emulsions, i.e. droplets (bubbles) of liquid (vapor) dispersed in a continuum liquid phase, typically water. In these multiphase systems, the transition between the dense and light phases is controlled by the surface tension which is, in turn, dictated by molecular interactions, notably the strength of the potential and its spatial range. That leads to interface widths of the order of the spatial range of such interactions, typically nanometers or below. A typical estimate of the interface width is as follows:
where denotes the surface tension. Typical values, in MKS units are and , deliver , i.e a fraction of nanometer.
As shown in Fig. 9, such nanometric interfaces challenge the lowKnudsen assumption which lays at the foundations of the hydrodynamic description. In fact, given that the interface thickness is comparable with the molecular mean free path, the result are local Knudsen numbers of the order unity. In the last decade, a number of technical extensions of the original LB method have been developed with the aim of gaining insights into these complex nonequilibrium interfacial phenomena (and yet, much remains to be done).
Here, we limit ourselves to summarize the main upgrades, namely: i) HigherOrder lattices, ii) Kinetic Boundary Conditions, iii) Regularization techniques.
iv.3.1 HigherOrder Lattices
As discussed in the initial part of this paper, hydrodynamics represents the “infrared” limit of kinetic theory, whereby all macroscopic heterogeneities live on much longer scales than the molecular ones (hydrodynamic transport regime). From the formal viewpoint, this means that the Boltzmann distribution is accurately described by its lowerorder moments, typically density (order ), current (order ) and momentumflux (order ).
In the hydrodynamic transport regime, all higher order moments (nonequilibrium excitations) are directly enslaved to the low order ones, hence they have no independent dynamics of their own.
Far from equilibrium, where strong inhomogeneities may persist down to nearmolecular scales, such loworder picture breaks down, and more kinetic moments concur to define the motion of matter beyond the hydrodynamic regime.
Without entering details, it is intuitively clear that this far from equilibrium regime requires the use of higher order lattices, securing the proper recovery of the correspondingly higher order moments. A formal theory of LB beyond NS, based upon higher order lattices, was laid down by X. Shan et al Shan et al. (2006). Those authors developed the discrete analogue of Grad’s expansion in Hermite polynomials, with full details on its specific implementation on a series of higher order lattices associated with different numerical quadrature rules. The first observation as compared to Grad’s 13 moment formulation is that the scheme provides a larger set of degrees of freedom Grad (1949). The discrete speeds of the corresponding lattices are typically of order or more, hence many more than the Grad’s moments Meng and Zhang (2011); Zhang et al. (2006).
iv.3.2 Kinetic Boundary Conditions
Higherorder lattices offer room for extramoments, but this is not sufficient perse, unless suitable boundary conditions are formulated near solid boundaries. This is where the lattice formulation makes a distinct contribution: while it appears very hard to devise wellposed boundary conditions for the kinetic moments, in fact complicated nonlinear tensors, lattice formulations lend themselves to conceptually transparent and mathematically wellposed formulations. The reason is always the same: the information moves along straight lines and boundary conditions can be formulated in terms of mechanical relations between the outgoing (fluidtowall) and the incoming (walltofluid) discrete distributions. Nonequilibrium flows exchange momentum with the solid walls along both tangential and normal directions: describing this exchange is the mandate of Kinetic Boundary Conditions. Empirical forms adapting this constraint to the lattice environment have been formulated in terms of accommodation coefficients Succi (2002b). Following Maxwell, the idea is that molecules impinging on the wall loose track of their incoming speed Maxwell (1878). Consequently, they are reinjected into the fluid along a random direction, with a velocity drawn from a Maxwellian at the local wall speed and temperature. Albeit handy, these accommodation schemes remain empirical in nature. A more satisfactory formulation was developed by Ansumali and Karlin, who basically expressed the accommodation coefficients in terms of the outgoing (fluidtowall) distribution functions and their equilibrium version, thus providing a closed and consistent recipe Ansumali and Karlin (2002). The AnsumaliKarlin boundary conditions exhibit a number of appealing properties. First, they preserve the positivity of the distribution at boundary nodes, i.e., if the incident distribution is positive, the reflected one is guaranteed to be positive too. Second, they readily extend to general wall scattering kernels, such as those allowing a blend of slip and reflection. At the time of this writing, these kinetic boundary conditions are the tool of the trade for finiteKnudsen LB simulations.
iv.3.3 Regularization
Regularization is a general and powerful idea across many fields of theoretical physics, to remove various forms of divergences and singularities which arise whenever a given description/theory fails to capture the physics in point.
Kinetic theory is no exception. Indeed, it is known since long that posthydrodynamic equations beyond the NS level, suffer a number of problems, primarily shortscale linear instabilities~Bobylev (1982). Many regularization schemes have been proposed ever since to tame such instabilities Struchtrup and Torrilhon (2003, 2007). Regularization procedures have been (re)discovered only recently by the LB community Latt and Chopard (2006); Montessori et al. (2014) and it is not yet clear how they relate to the corresponding counterparts in continuum kinetic theory. LB regularization consists of filtering out the contribution of nonhydrodynamic moments (ghosts) to the hydrodynamic ones: mass, momentum and momentumflux.
Let us dig a little bit deeper into the subject.
The standard LB scheme in BGK form reads (timestep made unit for simplicity):
(31) 
The actual distribution can be split as follows:
where the hydrocomponent collects terms up to third order Hermite polynomials associated with mass, momentum, momentumflux and energy flux, while collects all higher order terms, transport plus genuinely kinetic fields with no immediate macroscopic interpretation (ghosts in LB jargon). Formally, one defines a regularisation operator , projecting the actual distribution onto the hydrodynamic subspace, i.e., , that is .
Both hydro and ghost terms further split into equilibrium and nonequilibrium components, the ghost equilibrium being identically null by construction. Thus, what remains to be filtered out is just the ghost nonequilibrium, which is constantly revived at each freestreaming step, the nonequilibrium engine.
By applying the regularisation projector to the righthandside of the BGK equation, an operation corresponding to filtering out ghost components after streaming, we obtain the following RegularizedLB:
(32) 
The RegularizedLB has recently made proof of providing significant benefits in terms of improving the stability of the LB scheme under finiteKnudsen conditions.
The “posthydro” LB literature is vast and growing, but not conclusive yet. In particular, it is not clear how the three main ingredients mentioned above should be combined in order to obtain correct finiteKnudsen behavior.
A very valuable step in the comprehension of LB posthydro capabilities was provided by Ansumali, Karlin and coworkers, who discovered exact solutions to the hierarchy of nonlinear LB kinetic equations for stationary planar Couette flow at nonvanishing Knudsen numbers Ansumali et al. (2007). By using a 16speed twodimensional lattice and kinetic boundary conditions, these authors have derived closedform solutions for all higherorder moments, and solved them analytically.
The results indicate that the LB hierarchy with larger velocity sets does indeed approximate kinetic theory beyond the NS level. If only for a simple setup, those exact solutions indicate that LB equipped with kinetic boundary conditions is able to carry quantitative nonhydrodynamic information, hence it can be regarded as a kinetic closure in its own right.
It thus appears that the extension of the LB formalism to higher order lattices can result into an effective tool to probe deeper into nonequilibrium regimes beyond the hydrodynamic description Gan et al. (2015).
iv.4 Fluctuating Lattice Boltzmann
For the case of nanoscale flows, reproducing thermal behavior implies that the LB must incorporate the effects of statistical fluctuations.
To achieve this goal, following in the footsteps of LandauLifshitz fluctuating hydrodynamics, Ladd added a source of random fluctuations of the momentumflux tensor Ladd (1993), namely:
(33) 
where is the shear tensor and is the amplitude of the fluctuations.
The latter must be tuned so as to comply with the FluctuationDissipationTheorem. The resulting fluctuating NS equation reads
(34) 
where is the fluctationfree momentumflux tensor and is the body force acting on the fluid.
The LB scheme is modified accordingly, by adding a stochastic source to the right hand side:
(35) 
as sketched in Fig. 10. The source term is local in space and time and acts at the level of the stress tensor and nonhydrodynamic modes, with no effect on mass and momentum conservation.
cd
In actual practice, the source term is constructed via a set of the lattice eigenvectors with orthonormal according to the scalar product . In the D3Q19 scheme, the eigenvectors correspond to the kinetic moments: is relative to density, to current, to the momentum flux tensor, the remaining eigenvectors to nonhydrodynamic modes. The stochastic forcing reads as follows:
(36) 
where is a set of random numbers with zero mean and unit variance. Given the fact that the thermal velocity is fixed and equals the underlying lattice speed , the thermal mass is chosen in such a way as to obtain the thermal fluctuations according to .
The stochastic forcing has been subsequently improved so as to produce consistent fluctuations at all spatial scales, in particular at short distances where the effect on the translocating molecule is critical Dünweg et al. (2007); Adhikari et al. (2005). The fluctuating LB passes a number of litmus tests, particularly the compliance of velocityvelocity and forceforce autocorrelation functions with the principle of stochastic kinetic theory, in particular the fluctuationdissipation theorem.
iv.5 Lattice Boltzmann for nonideal fluids
Boltzmann originally derived his equation under the assumption of diluteness, whereby molecular collisions take place as zeroranged, instantaneous events, leading to large scattering angles. This rules out longrange, soft interactions giving rise to smallangle deflections. Such interactions, however, are of utmost importance for biological applications, since soft interlayers based on membranes and biopolymers define the spatial boundaries between different phases in biological systems. In aqueous media containing various ions, interactions are governed by a complex interplay of generic and specific interfacial interactions, typically controlled by dispersion and electrostatic forces, and locally shaped as hydrogen bond networks, disulphide bridges, hydrophobic, entropyinduced interactions and so on.
Softcore interactions can be included in the kinetic equation in the form of effective onebody forces, resulting from the collective interaction of a representative particle with the selfconsistent environment (bath) due to all other particles. Formally, such effective onebody force takes the following form
(37) 
where is the (unknown) twobody distribution and is the twobody atomistic potential. Following a customary practice, one writes , where , the twobody correlation function, collects the twobody physics. The exact form of the correlation function is known exactly only in a few precious cases; yet one can introduce several useful ansatzs which turn the formal expression 37 into an operational one. A typical ansatz, which has proven very fruitful for LB modelling of nonideal fluids, looks like follows:
(38) 
where is a model Green function and is a local functional of the fluid density . The above expression is quite general and allows a wide degree of latitude in modeling nonideal fluid interactions.
iv.6 Pseudopotential models
The most popular lattice transcription of the above, due to Shan and Chen Shan and Chen (1993, 1994), reads as follows:
(39) 
The sum runs over the prescribed set of neighbors, typically the first Brillouin region in the original version, and subsequently extended to the second or even the third one. Typically, all discrete speeds in the same Brillouin region share the same , so that in the original ShanChen (SC) formulation, there is just one coupling strength (see Fig.11). This is nonetheless sufficient to generate the main ingredients of nonideal fluids, namely a nonideal equations of state supporting phasetransitions, as well as surface tension. Moreover, this limitation can be readily lifted by extending the ShanChen interaction beyond the first Brillouin region, thereby explicitly accounting for both repulsive and attractive interactions.
The ShanChen expression delivers a nonideal equation of state of the form
(40) 
By choosing the generalized density in the form:
it is readily checked that the ShanChen fluid becomes critical for at a critical density . Note that even though the interaction is purely attractive (), no pileup instabilities take place because the force becomes increasingly faint as the density increases. This is the reason for using the generalized density instead of the physical one . This is a very expedient trick to trigger phasetransitions without any repulsive force. The translation from the force to the LB source proceeds through a systematic expansion in lattice Hermite polynomials. To leading order , but higher order terms, at least quadrupole ones, are definitely needed to obtain accurate results.
More generally, thick interfaces lead to large values of the Cahn number, namely the ratio of the interface width to a typical mesoscopic scale, say the droplet (bubble) diameter:
For most applications is in the range of tens, up to one hundred, microns, leading to very small Cahn numbers, . Replicating this number in LB simulations is totally unviable, for it would amount to placing lattice spacings across the particle diameter. LB simulations are forced to operate at much higher Cahn numbers of the order of , which means that potential inaccuracies due to finite Cahn number effects need to be carefully monitored.
Lack of sufficient symmetry also leads to the appearance of spurious currents, which may seriously affect the physical results whenever the density ratio between the light and dense phases is above 30 – 50, a problem which is further exacerbated in the case of fastmoving interfaces. Finally, the model is not thermodynamically consistent, as it is not derived from a meanfield freeenergy functional. That limitation is, however, less serious than it seems, since the method is, actually, equipped with a quasifree energy Sbragaglia et al. (2009). As a matter of fact, more disturbing, instead, is the fact that the ShanChen state equation features a sound speed smaller in the liquid with respect to the vapour phase, with significant consequences on the stability of the light phase across the interface.
Most of these limitations have been significantly mitigated, if perhaps not fully resolved, by subsequent developments. Among others, particularly noteworthy for the simulation of softglassy materials, is the socalled multirange pseudo potential method.
The main idea is to augment the original ShanChen shortrange attraction in the first Brillouin region (belt, in LB jargon), say the D3Q27 lattice in three dimensions, with a repulsive interaction acting on the second Brillouin region, i.e., including discrete velocities up to . The main advantage of this formulation is that a proper tuning between the firstbelt attraction and the secondbelt repulsion, makes possible to achieve smaller values of the surface tension, thereby permitting to sustain longlived, multidroplet configurations with a highly complex interfacial dynamics. Despite its many limitations, the SC method remains the most popular version of LB for nonideal fluids, mostly on account of its simplicity and transparency.
However, for applications leading to complex interfacial dynamics, more advanced schemes are required.
iv.6.1 FreeEnergy models
Another successful route to LB schemes for nonideal fluids is (lattice) Density Functional Theory (DFT), whereby the nonideal physics is encoded within a freeenergy functional of the form:
(41) 
where is a local functional of the fluid density and its gradients Swift et al. (1996).
The former encodes the nonideal equation of state, while the latter describes the effect of surface tension.
Variational minimization over density changes delivers the equations of motion of the nonideal fluids, governed by the Korteweg pressure tensor:
(42) 
where
is the nonlocal pressure, consisting of the bulk contribution , fixing the equation of state, plus an interface contribution associated to surface tension. In the above, stands for the space derivative along direction and is a tunable coefficient controlling the surface tension.
From the practical standpoint, the freeenergy formulation leads to nonlocal equilibria, which involve second order derivatives of the density field, thus taxing the simplicity of the method, and sometimes its stability as well, as compared to the pseudopotential method. Nevertheless, the freeenergy method has found broad use for many applications such as microflows over geometrically or chemically patterned surfaces and various types of droplet motions Liu et al. (2014).
Several variants have been developed over the years, which have considerably improved over the original versions, especially in terms of reaching higher density ratios without compromising numerical stability. Among others, highorder finitedifference schemes Lee and Lin (2005) and mergers with fluxlimiting methods Sofonea et al. (2004). These methods are currently used to simulate a variety of multiphase and multicomponent fluid applications, such as RayleighTaylor instabilities, droplet collisions, cavitation and freesurface flows. For a review, see the recent books Succi (2018); Montessori and Falcucci (2018) and references therein.
iv.6.2 Chromodynamic models
Another class of LB methods for nonideal fluids which has been revamped in the recent past is the socalled ColorGradient model. Here, the two phases/components, call them Red and Blue for convenience, segregate based on a topdown “colouring” rule, which sends the particles along the color gradient, namely Red towards Red and Blue towards Blue, thereby triggering an instability leading to interface formation Gunstensen et al. (1991).
The recoloring amounts to correcting the the Red and Blue populations as follows:
(43) 
where is the mass density fraction, is the cosine of the angle between the color gradient and the th discrete velocity. In the above , where star indicates the population after the application of the nonideal force stemming from surface tension, and , where superscript denotes the equilibria at zero flow. Finally, is a free parameter controlling the strength of the recoloring procedure, to all effects and purposes an antidiffusive operator. The crucial term is the second one at the righthand side of eq. (43), which, by construction, is active only at the interface between the two components. For more details see Leclaire et al. (2017) and Montessori et al. (2018a).
The interface is then stabilised by means of a chromodynamic force proportional to the amplitude of the color gradient, but opposite to it, so as to level out the deficit of one species over the other (color gradient). Although essentially ruledriven, modern variants of this scheme have proved capable of accessing parameter regimes which appear to be offlimits for ShanChen schemes and extensions thereof, as well as for FreeEnergy methods.} For instance, such methods have been recently applied to the design of microfluidic devices for droplet generation, such as flowfocusers and stepemulsifiers \cite{AM2018a,AM2018b}.
The interface is then stabilised by means of a chromodynamic force proportional to the amplitude of the color gradient, but opposite to it, so as to level out the deficit of one species over the other (color gradient). Although essentially ruledriven, modern variants of this scheme have proved capable of accessing parameter regimes which appear to be offlimits for ShanChen schemes and extensions thereof, as well as for FreeEnergy methods.
The recoloring amounts to correcting the Red and Blue populations as follows: where is the mass density fraction, is the cosine of the angle between the color gradient and the th discrete velocity. In the above , where the star indicates the population after the application of the nonideal force stemming from surface tension, and , where superscript denotes the equilibria at zero flow. Finally, is a free parameter controlling the strength of the recoloring procedure, to all effects and purposes an antidiffusive operator. The crucial term is the second one at the righthand side of eq. (43), which, by construction, is active only at the interface between the two components. For more details see Leclaire et al. (2017); Montessori et al. (2018a).
The interface is then stabilised by means of a chromodynamic force proportional to the amplitude of the color gradient, but opposite to it, so as to level out the deficit of one species over the other (color gradient). Although essentially ruledriven, modern variants of this scheme have proved capable of accessing parameter regimes which appear to be offlimits for ShanChen schemes and extensions thereof, as well as for FreeEnergy methods. For instance, such methods have been recently applied to the design of microfluidic devices for droplet generation, such as flowfocusers and stepemulsifiers Montessori et al. (2018b, c).
iv.6.3 Entropic models for multiphase flows
The lattice Boltzmann method can also be formulated by minimizing a suitable lattice Hfunction of the form (KullbackLeibler entropy):
where are suitable lattice weights. The resulting scheme takes the usual form of a standard LB, with the crucial twist that the relaxation time is adaptively adjusted in such a way as to secure compliance with the second principle of thermodyamics, namely:
Leaving a detailed description to the original literature, here we just mention that compliance with the second principle translates into a significant enhancement of numerical stability Karlin et al. (1998). For this reason, ELB has found profitable use for the simulation of lowviscous flows typical of fluid turbulence.
In a recent time, the ELB has been extended to the case of multiphase flows, and shown to provide stability benefits also in the viscous regime of relevance to many biological applications Wöhrwag et al. (2018); Montessori et al. (2017). Although these developments are too recent to permit a solid prononciation, the results appear very encouraging and leave hope that the entropic version of LB may become a major player in the field for the years to come.
V Fluids and Particles: The Lattice Boltzmann  Particle Dynamics Scheme
The multiphase LB schemes discussed in the previous Section have generated a mainstream of applications in soft matter research, since they permit to deal with flows of major dynamic and morphological complexity, such as foams and emulsions. However, they are unsuited to handle rigid bodies suspended in the continuum phase, nor can they describe in detail mechanical properties of deformable objects, say membranes, vesicles, cells and other biological bodies. To this purpose, LB needs to be explicitly coupled to particle methods, tracking the dynamics of the biological bodies immersed in the flow.
Indeed, most flows of biological interest consist of biological bodies of assorted nature: cells, polymers, proteins, floating in a fluid solvent, typically water. Such flows often operate at low, often nearzero, Reynolds number Purcell (1977), but this does not mean that hydrodynamic interactions (HI) can be neglected. To the contrary, HI have repeatedly been shown to accelerate a variety of nanoscale biological transport processes, such as biopolymer translocation across nanopores, amyloid aggregation of proteins in the cell and other related phenomena. This explains why the combination of LB with particle dynamics has made the object of extensive research and the development of multiple simulation schemes.
v.1 The extended particle model (EPM)
Flows with suspended objects have been modelled since the early days of LB research, starting with the trailblazing work of A. Ladd Ladd (1994a, b). Ladd’s original method consists of tracking the motion of rigid spherical bodies under the impact of the surrounding solvent hitting the surface of the body. The scales pertaining to body and solvent are fairly separate, since the former is much larger and heavier than the latter, hence the mass ratio serves as a suitable scale separator.
In the EPM, the exchange of momentum between particles and LB fluid is a boundarycollision method whereby the suspended body interact with each other only through the mediation of local collisions with the solvent molecules (LB). The solventbody collision is conservative, the momentum lost (gained) by the solvent to the body is gained (lost) by the body on the solvent. The method therefore is based on the local exchange of momentum by computing the amount of momentum that every population hitting the body surface exchanges with the latter. Being the total force and torque the sum of fluidbody momentum exchange (or drag interactions), direct particleparticle mechanical interactions or, in case of microscopic conditions, stochastic forces due to the Brownian motion, the body linear momentum and angular momentum obey the classical equations of motion
(44) 
It is worth mentioning that the fluidparticle coupling is hydrokinetic in nature because it is treated at collisional level rather than at hydrodynamic level, so that, hydrodynamic forces such as longtime tails, naturally emerge at times larger than the LB marching one. The approach takes naturally into account the solvent fluctuations, if present, since the latter are transmitted across the body surface analogously to the drag forces (see Fig.12).
Since the solventbody collisions are conservative, no extra stochastic source is needed on the body side.
Ladd’s method finds its best use in simulating colloidal suspensions of rigid particles and in this sense has found a comparatively limited application in the biological context. However, it provides a first and reliable example of coupling between LB and suspended bodies. At the same time, the method shows some inaccuracy due to the lattice discreteness, especially for nearcontact colloidcolloid interactions, occurring below the grid scale where the LB method can not resolve the lubrication forces. Some variants have been developed in the literature, making use of grid refinement and/or dynamic interpolations. Unsurprisingly, they only add to the computational complexity of the method, which is comparatively laborious on its own, as it demands full knowledge of the local fluidsolid connectivity, namely the dynamic list of fluid nodes interacting with the solid ones.
v.2 The pointparticle model (PPM)
A qualitatively different strategy has been proposed by Ahlrichs and Dünweg Ahlrichs and Dünweg (1998). This is a minimal way to embed pointlike particles in the LB fluid which, as opposed to Ladd’s method, is entirely forcedriven, hence metric instead of topologic. In essence, particles carry a phenomenological friction coefficient (also known as the bare friction coefficient) and the fluid exerts a drag force based on the instantaneous difference between the particle and fluid velocities, reading
(45) 
The addition of the point force into the fluid equations introduces a singularity into the flow field, as represented by Dirac’s delta function. This is automatically smoothed out by interpolation procedures, but requires nonetheless due care in the numerical treatment.
On the other hand, the flow field around a finitesized particle is generated by a distributed force located around the particle position, as shown in Fig. 13.
In the scheme, the flow field is treated as everywhere finite and the force density acts onto the LB fluid at position as
(46) 
where is a function interpolating between the particle position and the surrounding mesh node . In its original formulation, the method replaced and by the values of these fields at the mesh node nearest to the particle position.
A refined version resorts to a simple linear interpolation scheme using the mesh points on the elementary lattice cell containing the particle. Denoting the relative position of the particle in this cell by with the origin being at the lower left front edge, one defines , etc., the interpolation weights. The interpolated velocity then reads , where the sum is over the mesh points on the considered elementary lattice cell. The force exerted back to the fluid can then be chosen with the same weight coefficients or by spreading the force equally on the edges of the cell, the corresponding rule to preserve linear momentum being easily derived.
Due to the regularization of the pointforce, one should expect that the mobility of the suspended particle is not simply given by the bare friction coefficient , but is somehow renormalized. Indeed, the mobility is given by the sum of the bare mobility and a Stokestype contribution due to the lattice discretization. The effective friction coefficient relates to the bare one via , where is a geometrical correction coefficient Ahlrichs and Dünweg (1998).
The PPM finds application in simulating microscopic systems, with a stochastic term added on each particle in addition to the fluctuating LB bath. In this way, the particle does not leak energy on average and the fluctuationdissipation balance maintains a welldefined temperature of the fluidparticle system. Again, to balance and preserve momentum, such stochastic force is also restituted to the fluid with an opposite sign. The presence of longtime tails, that is, the inherent hydrodynamicallysustained motion of the moving particle generating a longtime decay of velocity, has also been observed Ahlrichs and Dünweg (1998). On the downside, the PPM is permeable to fluid momentum, as the local force is unable to create enough resistance to the incoming flow, and the classical Stokeslike picture of the streamlines does not apply.
Owing to its inherent simplicity, the PPM was first recognized to be useful to simulate topologically connected particles, such as polymers, in the presence of hydrodynamic interations. Whenever the embedded particles represent micro/mesoscopic objects, such as atoms or the beads of a polymer, mass diffusion plays a role comparable to mechanical and hydrodynamic forces. Another interesting extension involves coupling PPM with the Shan–Chen multicomponent LB to simulate efficiently complex fluid–fluid interfaces. The idea is to introduce a solvation freeenergy for the particle–fluid interaction proportional to the fluid density gradients Sega et al. (2013) and reads as follows
(47) 
Such force drives particles towards maximal or minimal density gradients, depending on the sign of the coupling coefficient . The approach is particularly appealing when used in conjunction with multiphase conditions. As a matter of fact, it makes possible to treat different multiphase fluids in the presence of suspended molecules, such as amphiphilic molecules as surfactants or, by adding another level of detail with electrostatics, polyelectrolytes in bicomponent fluids.
v.3 The diffusedparticle model (DPM)
An extension of PPM represents particles with finitesize extension, still relying on a forcebased mechanism. Given the finite extension, it is well suited for anisotropic particles, providing a very handy computational flexibility for the description of biological suspensions, cellular compartments or even entire cells (recalling that the interior of the cell is anisotropic due to intracellular organelles). The strategy is to consider the particle rototranslational response as originating from the coupling of the finitesize extension of the particle with the fluid momentum and vorticity. The particle is an effective diffused body, with no need to track its boundary to control its coupling with the environment.
The particle shape is described as an ellipsoid having three major radii along the three principal directions, with . To fit in the discrete nature of the lattice, are taken as three integers, one for each cartesian component (such requirement can be lifted but a few interesting properties described below, would be lost). The rototranslation is governed by rigid body dynamics of eqs. 44. The particle rotational state is encoded by the matrix whose rows are three orthogonal unit vectors aligned along the principle axis of the particle, that is, the basis to transform between the laboratory and the moving frames. The rototranslational state is given by the tensorial product
(48) 
where
is a shape function having compact support and the normalization property . Typically corresponds to a spherically symmetric diffused particle with a support extending over 64 mesh points. The translational response of the suspended body is designed according to the fluidparticle exchange kernel
where and the hydrodynamic DPM force is obtained via integration over the particle spatial extension. It reads as follows:
where .
The coupling between the body motion and the fluid vorticity is represented by a rotational kernel
with a rotational coupling coefficient. The corresponding drag torque is
where