A unified operator splitting approach for multi-scale fluid-particle coupling in the lattice Boltzmann method
Ulf D. Schiller
% A unified framework to derive discrete time-marching schemes for coupling of immersed solid and elastic objects to the lattice Boltzmann method is presented. Based on operator splitting for the discrete Boltzmann equation, second-order time-accurate schemes for the immersed boundary method, viscous force coupling and external boundary force are derived. Furthermore, a modified formulation of the external boundary force is introduced that leads to a more accurate no-slip boundary condition. The derivation also reveals that the coupling methods can be cast into a unified form, and that the immersed boundary method can be interpreted as the limit of force coupling for vanishing particle mass. In practice, the ratio between fluid and particle mass determines the strength of the force transfer in the coupling. The integration schemes formally improve the accuracy of first-order algorithms that are commonly employed when coupling immersed objects to a lattice Boltzmann fluid. It is anticipated that they will also lead to superior long-time stability in simulations of complex fluids with multiple scales. %The implementation of the coupling schemes is straightforward, and it %is anticipated that long-time stability in simulations of complex %fluids will be improved. %
lb-operator-splitting.tex
\documentclass[preprint,
               sort&compress,
              ]{elsarticle}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{mathtools}
\usepackage{algorithm}
\usepackage{algorithmicx}
\usepackage{algpseudocode}
\usepackage{hyperref}
\usepackage[normalem]{ulem}
\usepackage{xcolor}

%\newcommand{\comment}[1]{{\color{red}(#1)}}
%\newcommand{\new}[1]{{\color{blue} #1}}
%\newcommand{\cancel}[1]{{\color{green}{{#1}}}}
%\newcommand{\replace}[2]{{\color{green}{{#1}}} $\Rightarrow$ {\color{blue}#2}}
%
%\renewcommand{\comment}[1]{{}}
%\renewcommand{\new}[1]{#1}
%\renewcommand{\cancel}[1]{{}}
%\renewcommand{\replace}[2]{#2}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\renewcommand{\vec}[1]{{\ensuremath{\mathchoice
                     {\mbox{\boldmath$\displaystyle\mathbf{#1}$}}
                     {\mbox{\boldmath$\textstyle\mathbf{#1}$}}
                     {\mbox{\boldmath$\scriptstyle\mathbf{#1}$}}
                     {\mbox{\boldmath$\scriptscriptstyle\mathbf{#1}$}}}}}%

\newcommand{\tens}[1]{{\ensuremath{\mathchoice
                     {\mbox{$\displaystyle\mathsf{#1}$}}
                     {\mbox{$\textstyle\mathsf{#1}$}}
                     {\mbox{$\scriptstyle\mathsf{#1}$}}
                     {\mbox{$\scriptscriptstyle\mathsf{#1}$}}}}}%

\renewcommand{\d}{\ensuremath{\partial}}

\newcommand{\Dt}{\ensuremath{\Delta t}}
\newcommand{\eq}{\ensuremath{\text{eq}}}
\newcommand{\nq}{\ensuremath{\text{neq}}}
\newcommand{\Fvec}{\ensuremath{\vec{F}}}
\newcommand{\Rvec}{\ensuremath{\vec{R}}}
\newcommand{\Uvec}{\ensuremath{\vec{U}}}
\newcommand{\Vvec}{\ensuremath{\vec{V}}}
\newcommand{\Xvec}{\ensuremath{\vec{X}}}
\newcommand{\avec}{\ensuremath{\vec{a}}}
\newcommand{\cvec}{\ensuremath{\vec{c}}}
\newcommand{\fvec}{\ensuremath{\vec{f}}}
\newcommand{\pvec}{\ensuremath{\vec{p}}}
\newcommand{\rvec}{\ensuremath{\vec{r}}}
\newcommand{\uvec}{\ensuremath{\vec{u}}}
\newcommand{\vvec}{\ensuremath{\vec{v}}}
\newcommand{\xvec}{\ensuremath{\vec{x}}}
\newcommand{\zetavec}{\ensuremath{\vec{\zeta}}}
\newcommand{\Pitens}{\ensuremath{\tens{\Pi}}}

\newcommand{\1}{\ensuremath{\mathsf{I}}}
\newcommand{\Iop}{\ensuremath{\mathcal{I}}}
%\newcommand{\Lop}{\ensuremath{\mathcal{L}}}
\newcommand{\Lambdaop}{\ensuremath{\mathsf{\Lambda}}}
\newcommand{\Omegaop}{\ensuremath{\mathsf{\Omega}}}

\def\bibsection{\section*{References}}

%\algrenewcommand\Require{\algorithmicrequire\ }%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{A unified operator splitting approach for multi-scale
fluid-particle coupling in the lattice Boltzmann method}

\author{Ulf D. Schiller}\fnref{fn1}

\ead{u.schiller@ucl.ac.uk}

\address{%
  Theoretical Soft Matter and Biophysics,
  Institute of Complex Systems,
  Forschungszentrum~J\"ulich,
  52425 J\"ulich,
  Germany
}

\fntext[fn1]{Present address: Centre for Computational Science, University College London, 20 Gordon Street, London WC1H 0AJ, United Kingdom}

%\address{%
%  Department of Chemical Engineering,
%  University of Florida,
%  Gainesville, 
%  Florida 32611-6005,
%  USA
%}

\date{\today}



\begin{abstract}
%
A unified framework to derive discrete time-marching schemes for
coupling of immersed solid and elastic objects to the lattice
Boltzmann method is presented. Based on operator splitting for the
discrete Boltzmann equation, second-order time-accurate schemes for
the immersed boundary method, viscous force coupling and external
boundary force are derived. Furthermore, a modified formulation of the
external boundary force is introduced that leads to a more accurate
no-slip boundary condition. The derivation also reveals that the
coupling methods can be cast into a unified form, and that the
immersed boundary method can be interpreted as the limit of force
coupling for vanishing particle mass. In practice, the ratio between
fluid and particle mass determines the strength of the force transfer
in the coupling. The integration schemes formally improve the accuracy
of first-order algorithms that are commonly employed when coupling
immersed objects to a lattice Boltzmann fluid. It is anticipated that
they will also lead to superior long-time stability in simulations of
complex fluids with multiple scales.
%The implementation of the coupling schemes is straightforward, and it
%is anticipated that long-time stability in simulations of complex
%fluids will be improved.
%
\end{abstract}



\begin{keyword}
lattice Boltzmann
\sep fluid-particle coupling
\sep force coupling
\sep immersed boundary method
\sep external boundary force
\end{keyword}



\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

The lattice Boltzmann method (LBM) \cite{Succi-LB} is a popular
approach to simulate hydrodynamic phenomena and is widely used as an
alternative to continuum-based approaches to solving the Navier-Stokes
equations. The LBM is a mesoscopic model that bridges between kinetic
theory and macroscopic hydrodynamics, i.e., it reproduces the
conservation laws of the Boltzmann equation (in particular mass and
momentum conservation) on macroscopic time and length scales
\cite{He1997a,He1997b,Shan1998,Junk2005a}.
%
Since the LBM enables simple and accurate incorporation of geometric
boundaries, it is widely applied to fluid-structure interaction
problems \cite{Krafczyk2001,Geller2006,DeRosis2013}. 
%
%In macroscopic flows, however, it may be faced with numerical
%instabilities that arise in the low-viscosity regime
%\cite{Brownlee2007}.  Additional recipes such as the entropic LBM
%\cite{Karlin1998,Boghosian2001,Succi2002} have been developed to
%preserve stability and enable, e.g., simulations of turbulent flows
%\cite{Chikatamarla2013}.

In recent years, there has been growing interest in using the lattice
Boltzmann method to simulate momentum transport and hydrodynamic
interactions in complex fluids and soft matter systems
\cite{Duenweg2008,Aidun2010}.  Such systems concern objects that are
immersed in a continuum solvent, for example, colloidal suspensions
\cite{Ladd1994a,Ladd2001}, polymer solutions
\cite{Ahlrichs1999,Pham2009,Ladd2009}, or biological cells
\cite{Alexeev2005,Dupin2007,Wu2010}. These objects can be represented
by a set of ``coarse-grained'' particles which interact via effective
potentials and reproduce the correct dynamics on physically
interesting scales. Since these scales are typically large compared to
atomistic scales, the solvent can be treated as a continuum and the
hydrodynamic fields represent slow degrees of freedom of the
microscopic kinetics.
%
The standard lattice Boltzmann algorithm gives a second-order accurate
approximation of the macroscopic transport equations for the conserved
moments, however, for grid-scale Reynolds numbers above unity the
non-conserved moments may be subject to oscillations that can cause
nonlinear instabilities \cite{Brownlee2007}.
%Additional recipes such as the entropic LBM
%\cite{Karlin1998,Boghosian2001,Succi2002} have been developed to
%preserve stability and enable, e.g., simulations of turbulent flows
%\cite{Chikatamarla2013}.
%
Soft matter flows are usually in the creeping flow regime, i.e., they
are incompressible (low Mach number) and inertia can be neglected (low
Reynolds number). 
%
%\replace{%
%%
%The accuracy-stability dilemma regarding grid-scale
%Reynolds number above unity \cite{Brownlee2007} thus rarely hampers
%soft matter simulations. The lattice Boltzmann method is therefore
%well suited to simulate hydrodynamics in this domain.
%%
%}{%
%%
%The stability conditions are typically satisfied in this regime, and
%the standard lattice Boltzmann algorithm is therefore well suited to
%simulate hydrodynamics in the soft matter domain.
%%
%}

%In macroscopic flows, however, it may be faced with numerical
%instabilities that arise in the low-viscosity regime
%\cite{Brownlee2007}.  Additional recipes such as the entropic LBM
%\cite{Karlin1998,Boghosian2001,Succi2002} have been developed to
%preserve stability and enable, e.g., simulations of turbulent flows
%\cite{Chikatamarla2013}.
%
%The lattice Boltzmann method is well suited to simulate hydrodynamics
%in this regime, and various approaches have been proposed to simulate,
%e.g., multiphase flows
%\cite{Shan1993,Swift1996,Luo1998,Gonnella1999,Luo2000a} and
%multi-component systems
%\cite{Luo2003,Guo2005,Arcidiacono2006,Halliday2007,Li2007a}.
%
%, particulate suspensions \cite{} and polymer systems \cite{}.
%
%The lattice Boltzmann method is a popular approach to simulate
%momentum transport and hydrodynamic interactions in complex fluids and
%soft matter systems \cite{Succi-LB,Duenweg2008,Aidun2010}. Typical
%examples of such systems are colloidal suspensions
%\cite{Ladd1994a,Ladd2001}, polymer solutions
%\cite{Ahlrichs1999,Pham2009,Ladd2009}, and more recently biological
%cells \cite{Alexeev2005,Dupin2007,Wu2010}. 

The dynamics of complex fluids with immersed objects can be strongly
influenced by the coupling between solute and solvent degrees of
freedom. The motion of a solute particle creates a perturbation of the
solvent that spreads out to other particles. Hence, the particle
movements become correlated due to viscous momentum diffusion. The
idea of hybrid fluid-particle methods is to simulate the momentum
transport in an explicit solvent with the LBM while the solute objects
are coupled to the flow field by some sort of interaction force
\cite{Duenweg2008}.
%The common idea of these computational approaches is to simulate the
%fluid degrees of freedom by the LBM, which is coupled to the dymamics
%of the solute objects by some sort of interaction. The role of the
%interaction is to enforce appropriate boundary conditions at the
%fluid-object interface. 
If the immersed objects have internal degrees of freedom,
e.g. polymers, membranes or cells, the equations of motion for these
degrees of freedom have to be solved simultaneously with the fluid
dynamics. Various approaches have been proposed to solve this task,
such as the immersed boundary method (IBM)~\cite{Peskin2002}, force
coupling through viscous drag~\cite{Ahlrichs1999}, or external
boundary force (EBF) \cite{Wu2010}. They have in common that a set of
Lagrangian particles or marker points is added to the lattice
Boltzmann fluid, which require the simultaneous integration of
Newton's equations of motion, e.g., by molecular dynamics (MD)
techniques.

Integration schemes for MD can be faced with spurious effects due to
the time discretization, for example, systematic drift of energy or
temperature, artificial spatial correlations between particles
\cite{Serrano2006}, and numerical instabilities \cite{Schlick1998}.  A
number of improved algorithms has been introduced,
e.g. \cite{Tuckerman1992,Ricci2003,Serrano2006,Bussi2007}, which are
commonly based on a Trotter expansion of the Liouville propagator
\cite{Trotter1959,Kato1978,Suzuki1985}. Integration schemes derived by
splitting methods can improve accuracy and lead to symplectic
algorithms that are usually also more stable
\cite{Tuckerman1992,McLachlan2002,Bou-Rabee2013}.

In this communication, a general approach to derive time-marching
schemes to solve hybrid fluid-particle systems is presented. It is
based on the interpretation of the lattice Boltzmann equation as a
Strang splitting scheme \cite{Dellar2013} combined with common
operator splitting techniques used in molecular dynamics (MD) and
related methods
\cite{FrenkelSmit,Tuckerman1992,Skeel1997,Shardlow2003,
DeFabritiis2006,Bou-Rabee2013}. The operator splitting makes it
straightforward to derive explicit time-discretizations, in
particular, second-order time-accurate integration schemes for the
various coupling methods are proposed. Furthermore, a modified
external boundary force is proposed that ensures the no-slip boundary
condition after application of the force operator. The latter also
reveals that the IBM can be interpreted as the limit of force coupling
methods for vanishing particle mass, and that the force transfer
during coupling is controlled by the ratio between fluid and particle
mass. This notion may also be applicable to other mesoscopic methods
such as multi-particle collision dynamics
\cite{Malevanets1999,Gompper2009}.

The remainder of this work is organized as follows: In section
\ref{sec:Boltzmann} and \ref{sec:splitting}, the derivation of the
lattice Boltzmann equation and its interpretation as an operator
splitting scheme \cite{Dellar2013} are briefly reviewed. In section
\ref{sec:CSC-LB}, an alternative lattice Boltzmann algorithm in the
form of a collide-stream-collide scheme is outlined and analyzed. The
operator splitting is then extended to several coupling methods, i.e.,
the immersed boundary method (section \ref{sec:IBM}), force coupling
(section \ref{sec:force-coupling}), and external boundary force
(section \ref{sec:ebf}). Finally, multiple time-step schemes are
briefly discussed in section~\ref{sec:multi-timestep}.



\section{Discrete Boltzmann equation}
\label{sec:Boltzmann}

In the kinetic theory of gases, the dynamics and transport properties
of fluids are determined by the Boltzmann equation
\cite{Cercignani1988}
%The lattice Boltzmann equation is a fully discretized form of the
%Boltzmann equation with a linear collision operator
%\cite{McNamara1988,Higuera1989,Benzi1992,Chen1998}
%
\begin{equation}\label{eq:Boltzmann}
%\left( \d_t + \vvec \cdot \d_\xvec + \avec \cdot \d_\vvec \right)
%f(\xvec,\vvec,t) = - \Omegaop ( f - f^\eq ) .
\left(\frac{\d}{\d t} + \cvec\cdot\frac{\d}{\d\xvec} +
  \avec \cdot\frac{\d}{\d\cvec}\right) f(\xvec,\cvec,t) =
- \Omegaop ( f - f^\eq ) . 
%\left(\d f\right)_\text{coll} .
\end{equation}
%
Here, $f(\xvec,\cvec,t)$ is the one-particle distribution function at
position $\xvec$, velocity $\vvec$, and time $t$, and the collision
operator $\Omegaop$ on the right-hand side has been linearized around
the equilibrium distribution $f^\eq$ \cite{Bhatnagar1954,Gross1959}.
%
In order to discretize the Boltzmann equation \eqref{eq:Boltzmann},
the moment integrals of the distribution function are formally
evaluated by a Gauss-Hermite quadrature, where the abscissae of the
quadrature comprise a set of discrete velocities $\cvec_i$
\cite{Shan1998,Shan2006}. This procedure leads to the {discrete
Boltzmann equation}, which is also the basis of discrete velocity
models \cite{Broadwell1964a,Luo2000}
%The lattice structure emerges from a set of discrete velocities
%$\cvec_i$ that link discrete positions on a spatial grid. The discrete
%velocities can formally be obtained as abscissae in a Gauss-Hermite
%quadrature of the moment integrals of the distribution function
%\cite{Shan1998,Shan2006}. This first discretization step leads to the
%discrete Boltzmann equation which is also the basis of so-called
%discrete velocity machines \cite{Broadwell1964a,Luo2000}.
%
\begin{equation}\label{eq:discrete-Boltzmann}
%\left( \d_t + \cvec_i \cdot \d_\xvec \right) f_i = - \sum_j \Omega_{ij} \left( f_j - f_j^{\eq} \right) + G_i ,
\left( \frac{\d}{\d t} + \cvec_i \cdot \frac{\d}{\d \xvec} \right) f_i = - \sum_j \Omega_{ij} \left( f_j - f_j^{\eq} \right) + G_i ,
\end{equation}
%
%The collision operator on the right hand side has been linearized
%around the equilibrium distribution $f_i^{\eq}$ \cite{Higuera1989},
where $G_i$ is a body force corresponding to the acceleration term
$-\avec\cdot\nabla_\cvec f$. 
%The discrete velocities $\cvec_i$ link
%discrete positions and thus induce a spatial lattice. 
The standard approach to space-time discretization of the discrete
Boltzmann equation is integrating along the characteristic $(\xvec+h
\cvec_i, t+h)$ and approximating the integral of the collision term by
the trapezium rule \cite{He1997a,He1997b}. This leads to a non-linear
implicit equation for $f_i(\xvec+h \cvec_i, t+h)$ that can be turned
into an explicit scheme by an appropriate re-definition of variables
\cite{He1998}. This explicit scheme is the lattice Boltzmann equation
in the standard form
%\cite{McNamara1988,Higuera1989,Benzi1992,Chen1998}
\cite{Benzi1992,Chen1998,He1998a,Guo2002}
%
\begin{equation}\label{eq:LBE}
\tilde f_i(\xvec+h\cvec_i,t+h) = \tilde f_i(\xvec,t) - \sum_j
\Lambda_{ij} \left( \tilde f_j - f_j^{\eq} \right) + \sum_j (\delta_{ij} - \frac{1}{2} \Lambda_{ij}) h G_j(t) ,
\end{equation}
%
where the new variables are \cite{He1998,Dellar2003}
%
\begin{equation}\label{eq:He-trafo1}
\tilde f_i(t) = f_i(t) + \sum_j \frac{h}{2} \Omega_{ij} \left( f_j - f_j^{\eq} \right) - \frac{h}{2} G_i(t) .
\end{equation}
%
It is important to make the distinction between $f_i$ and the new
variables $\tilde f_i$ and to replace the collision matrix $h\Omegaop$
by the discrete collision matrix $\Lambdaop =
(1+\frac{h}{2}\Omegaop)^{-1}h\Omegaop$
\cite{Henon1987,Xu2003,Dellar2003}. Without these replacements,
Eq.~\eqref{eq:LBE} is only a first-order accurate discretization of
the Boltzmann equation \cite{Ohwada1998}. While $\Lambdaop$ enters in
the collision step and leads to the correct viscosity, the variable
transformation \eqref{eq:He-trafo1} is essential to evaluate the
macroscopic flow fields which have to be obtained as moments of the
back-transformed variables $f_i$:
%
%However, the new variables $\tilde f_i$ and the collision matrix
%$\Lambdaop$ in \eqref{eq:LBE} do not correspond exactly to their
%counterparts in the Boltzmann equation \eqref{eq:discrete-Boltzmann}
%\cite{Xu2003}.  This is related to the fact that Eq.~\eqref{eq:LBE} is
%only a first-order accurate discretization \cite{Ohwada1998}, and the
%back-transformation of the variables $\tilde f_i$ is required to
%achieve second-order accuracy. 
%
%The macroscopic flow fields are obtained as moments of the
%back-transformed variables $f_i$
%
\begin{align}
\rho &= \sum_i f_i ,
& \rho \uvec &= \sum_i f_i \cvec_i ,
& \Pitens &= \sum_i f_i \cvec_i \cvec_i .
%& (\Pitens -\Pitens^\eq) &= \sum_i (f_i - f_i^\eq) \cvec_i \cvec_i \\
%&&& = \frac{1}{\rho} \sum_i \tilde f_i \cvec_i + \frac{h}{2} \Fvec
%&&  = \left( 1 - \frac{\Lambda_{ij}}{2} \right) \sum_i (\tilde f_i - \tilde f_i^\eq) \cvec_i \cvec_i 
\end{align}
%
In the absence of forces, the back-transformation is only relevant for
evaluating the (non-equilibrium) stress and does not change the
mass and momentum densities.

Dellar \cite{Dellar2013} has recently pointed out that the lattice
Boltzmann equation can also be interpreted as a Strang splitting
scheme. In the following, operator splitting for the discrete
Boltzmann equation is revisited, which will serve as the basis to
derive accurate coupling schemes in later sections.  The operator
splitting approach will also provide a complementary view on the role
of Eq.~\eqref{eq:He-trafo1}.




\section{Operator splitting for the discrete Boltzmann equation}
\label{sec:splitting}

In the operator splitting interpretation put forward by
Dellar \cite{Dellar2013}, the lattice Boltzmann equation is split into
decoupled streaming, collision and forcing steps. The idea is to find
approximate solutions to the individual steps and combine them into a
time-marching scheme, the accuracy of which depends on the order of
the applied splitting and the accuracy of the individual approximate
solutions.
%
Splitting methods \cite{McLachlan2002} make use of decomposition
formulas for exponential operators where the exponential of a sum of
non-commuting operators is approximated by a symmetrized product of
exponentials of the individual operators
\cite{Trotter1959,Kato1978,Suzuki1985}. For the simplest case of two
operators, the resulting splitting scheme is equivalent to the
difference scheme introduced by Strang \cite{Strang1968}.  
%
%\cancel{More generally, the composition of flows of the decoupled
%  differential equations can be written as products of exponential
%  operators \cite{Groebner1960}, which can be expressed as a single
%  exponential by application of the Baker-Campbell-Hausdorff formula
%  \cite{Varadarajan1974,Hairer2006}. The truncated Taylor series of
%  the exponentials can be applied to sufficiently differentiable
%  functions which leads to formally accurate splitting schemes
%  \cite{Yoshida1990,Hairer2006}.  }
%
%\replace{Strictly speaking, the decomposition formulas are valid for linear
%operators only, but it can be generalized to nonlinear operators with
%continuous second derivatives \cite{Bobylev2001}.}
%
The decomposition formulas are valid for linear operators and can be
generalized to nonlinear operators by application of the
Baker-Campbell-Hausdorff formula
\cite{Groebner1960,Varadarajan1974,Hairer2006}. The truncated Taylor
series of the exponentials can be applied to sufficiently
differentiable functions and yield formally accurate splitting schemes
\cite{Yoshida1990,Hairer2006}.
%
Hence it is justified to use splitting schemes for the Boltzmann
kinetic equation \cite{Bobylev2001} although the collision operator
depends non-linearly on the flow velocity through the equilibrium
distribution. The remainder of this section closely follows the work
of Dellar \cite{Dellar2013}. However, instead of combining adjacent
parts of different time-steps to recover the standard lattice
Boltzmann algorithm up to a transformation of variables, we will
combine the sub-parts of one time-step into an alternative update
scheme in section~\ref{sec:CSC-LB}. This provides a complementary view
on the splitting scheme and clarifies the role of the transformed
variables as intermediate values.

The discrete Boltzmann equation \eqref{eq:discrete-Boltzmann} can be
written in terms of operators as
%
\begin{equation}\label{eq:split-lb}
%
\frac{\d}{\d t} f_i(\xvec,t) = \left( \mathcal{S} + \mathcal{C} + \mathcal{F} \right) f_i(\xvec,t) ,
\end{equation}
%
where the operators for streaming $\mathcal{S}$, collisions
$\mathcal{C}$, and forces $\mathcal{F}$ are
%
\begin{subequations}\label{eq:decoupled}
\begin{align}
%
\mathcal{S} f_i &=  - \cvec_i \cdot \frac{\d}{\d \xvec} f_i , \\
%
\mathcal{C} f_i &= - \sum_j \Omega_{ij} \left( f_j - f_j^{\eq} \right) , \\
%
\mathcal{F} f_i &= G_i .
%
\end{align}
\end{subequations}
%
The action of forces is usually incorporated into the collision
operator $\mathcal{C}$, but it can also be treated as a third
decoupled step \cite{Salmon1999a,Salmon1999b,Dellar2013}. Here, we
adopt the splitting using three decoupled parts for streaming,
collisions and forces. The motivation is that the force operator
couples the interactions with the lattice Boltzmann equation and may
be implemented more cleanly if treated separately.

The formal solutions of the decoupled equations \eqref{eq:decoupled}
are
%
\begin{subequations}
\begin{align}
f_i(\xvec,t+h) &= e^{h\mathcal{S}} f_i(\xvec,t) \approx \mathsf{S} f_i(\xvec,t) , \label{eq:streaming} \\
f_i(\xvec,t+h) &= e^{h\mathcal{C}} f_i(\xvec,t) \approx \mathsf{C} f_i(\xvec,t) , \label{eq:collision} \\
f_i(\xvec,t+h) &= e^{h\mathcal{F}} f_i(\xvec,t) \approx \mathsf{F} f_i(\xvec,t) ,
\label{eq:force}
\end{align}
\end{subequations}
%
where the three operators $\mathsf{S} \approx \exp[h\mathcal{S}]$,
$\mathsf{C} \approx \exp[h\mathcal{C}]$, and $\mathsf{F} \approx
\exp[h\mathcal{F}]$ are introduced as approximations to the separate
exponential operators. Note that the index $i$ and the coordinate
$\xvec$ in the streaming part \eqref{eq:streaming} are to be
understood as ranging over all possible values because streaming is a
non-local operation. In this sense, the operator $\exp[h\mathcal{S}]$
acts on the whole ensemble of distributions.
%
The discrete approximations for the decoupled steps are obtained
below.

The formal solution to \eqref{eq:split-lb} is given by
%
\begin{equation}\label{eq:LB-evolution}
f_i(\xvec,t+h) = e^{h(\mathcal{S}+\mathcal{C}+\mathcal{F})}
f_{i}(\xvec,t) ,
\end{equation}
%
where again the index $i$ and the coordinates $\xvec$ are ranging over
all values.
%
A second-order Trotter-Suzuki decomposition
\cite{Trotter1959,Suzuki1985} of the exponential operator leads to the
splitting scheme
%
\begin{equation}\label{eq:LB-splitting}
\begin{split}
e^ { h ( \mathcal{S} + \mathcal{C} + \mathcal{F} ) } &= e^ {\frac{h}{2} \mathcal{F} } e^ { \frac{h}{2} \mathcal{C} } e^{ h \mathcal{S} } e^{ \frac{h}{2} \mathcal{C} } e^{ \frac{h}{2} \mathcal{F} } + O(h^3) , \\
&= \mathsf{F}^{\frac{1}{2}} \mathsf{C}^{\frac{1}{2}} \mathsf{S} \mathsf{C}^{\frac{1}{2}} \mathsf{F}^{\frac{1}{2}} + O(h^3) ,
\end{split}
\end{equation}
%
where the operators $\mathsf{S}$, $\mathsf{C}^\frac12$, and
$\mathsf{F}^\frac12$ are the discrete approximations of the
exponential operators for a step $h$ and $h/2$, respectively.
%
The lattice Boltzmann algorithm requires that the streaming part
$\mathsf{S}$ must not be split in order to be compatible with the
lattice structure. For second-order accuracy, $\mathsf{F}^\frac{1}{2}$
and $\mathsf{C}^\frac{1}{2}$ have to be arranged symmetrically but
their order is in principle arbitrary. The choice to put
$\mathsf{F}^\frac{1}{2}$ to the outside is motivated by the coupling
algorithms discussed below. For small enough forces
$\mathsf{F}^\frac{1}{2}$ and $\mathsf{C}^\frac{1}{2}$ commute,
cf. section~\ref{sec:forces}.




\subsection{Streaming}

The streaming part \eqref{eq:streaming} can be written as a derivative along the
characteristic $(x+h'\cvec_i,t+h')$
%
\begin{equation}
%( \d_t + \cvec_i\cdot\d_\xvec ) f_i(\xvec,t) = \frac{d}{d h'} f_i(\xvec+h'\cvec_i,t+h') = 0 ,
\left( \frac{\d}{\d t} + \cvec_i\cdot\frac{\d}{\d \xvec} \right) f_i(\xvec,t) = \frac{d}{d h'} f_i(\xvec+h'\cvec_i,t+h') = 0 ,
\end{equation}
%
and can be integrated from $0$ to $h$ to give
%
\begin{equation}
\begin{gathered}
\int_0^h \frac{d}{dh'} f_i(\xvec+h'\cvec_i,t+h') dh' = f_i(\xvec+h\cvec_i,t+h) - f_i(\xvec,t) = 0 .
\end{gathered}
\end{equation}
%
We thus get the discrete streaming step
%
\begin{equation}
\begin{gathered}
f_i(\xvec+h\cvec_i,t+h) = f_i(\xvec,t) ,
\end{gathered}
\end{equation}
%
from which the spatial grid structure emerges, i.e., only coordinates
connected by the discrete vectors $h \cvec_i$ are admissible. 

%Thus, strictly speaking, the coordinate $\xvec$ on the left hand sides
%of \eqref{eq:streaming} and \eqref{eq:LB-evolution} has to be replaced
%by $\xvec + h \cvec_i$ for the equalities to hold for given $i$ and
%$\xvec$.



\subsection{Collisions}

%\cancel{%
%%
%The discretization of the collision part \eqref{eq:collision} is
%obtained using the $O(h^3)$ accurate Crank-Nicolson rule
%\mbox{\cite{Dellar2003,Dellar2013}}
%%
%\begin{equation}\label{eq:CN-coll}
%\begin{split}
%\fvec(\xvec,t+h)-\fvec(\xvec,t) &= - \frac{h}{2} \Omegaop \left[ \fvec(\xvec,t+h) - \fvec^\eq(\xvec,t+h) + \fvec(\xvec,t) - \fvec^\eq(\xvec,t) \right] , \\
%\end{split}
%\end{equation}
%%
%where we have introduced the vector $\fvec=(f_0,f_1,\dots,f_q)^T$.
%%
%However, as was pointed out by Brownlee et al. \cite{Brownlee2007},
%$\fvec - \fvec^\eq$ can be of order $O(\tau^2)$ where $\tau$ is the
%collision time. Therefore, the accuracy of \eqref{eq:CN-coll} is in
%fact $O((h/\tau)^3)$, which can cause stability problems if the grid
%scale Reynolds number $h/\tau$ is above unity.
%%
%}

The collision part \eqref{eq:collision} leaves the equilibrium
distribution invariant since $\mathcal{C} f_i^\eq = 0$. Hence we
consider only the non-equilibrium part $f_i^\nq = f_i - f_i^\eq$ and
the differential equation
%
\begin{equation}
\frac{d}{d t} \fvec^\nq = \mathcal{C} \fvec^\nq = - \Omegaop \, \fvec^\nq ,
\end{equation}
%
where we have introduced the vector notation
$\fvec=(f_0,f_1,\dots,f_q)^T$.
%
An approximate solution is obtained using the $O(h^3)$ accurate
Crank-Nicolson rule \mbox{\cite{Dellar2003,Dellar2013}}
%
\begin{equation}\label{eq:CN-approx}
\begin{split}
  \fvec^\nq(t+h)-\fvec^\nq(t) &= - \frac{h}{2} \Omegaop \left[ \fvec^\nq(t+h) + \fvec^\nq(t) \right] , \\
\end{split}
\end{equation}
%
which leads to
%
\begin{equation}\label{eq:CN-coll}
\begin{split}
  \fvec^\nq(t+h) = \left( \1 + \frac{h}{2} \Omegaop \right)^{-1} 
  \left( \1 - \frac{h}{2} \Omegaop \right) \fvec^\nq(t) .
\end{split}
\end{equation}
%
A comparison with the exact solution 
%
\begin{equation}\label{eq:exact-coll}
  \fvec^\nq(t+h) = \exp [ - h \Omegaop ] \, \fvec^\nq(t)
\end{equation}
%
shows that the accuracy of \eqref{eq:CN-coll} is in fact
$O((h/\tau)^3)$, which can cause nonlinear instabilities for grid
scale Reynolds number $h/\tau$ above unity \cite{Brownlee2007}.  On
the other hand, the analysis of Dellar \cite{Dellar2013} revealed that
the exact solution \eqref{eq:exact-coll} leads to excessively rapid
decay of shear modes as $\tau \lesssim h$, while the Crank-Nicolson
approximation \eqref{eq:CN-coll} gives accurate decay rates even for
small collision times. The reason for this somewhat counter-intuitive
result lies in the compensation of errors introduced by the splitting
of the Boltzmann equation into uncoupled steps on one hand, and the
Crank-Nicolson approximation of the collision operator on the other.

%\cancel{
%%
%Using that $\fvec^{\eq}$ is invariant under collisions we can solve
%for $\fvec(\xvec,t+h)$
%%
%\begin{equation}\label{eq:LB-collisions}
%\begin{split}
%\fvec(\xvec,t+h) %&= \mathsf{C} \fvec(\xvec,t)\\
%%&= \fvec(\xvec,t) - \left( \1 + \frac{h}{2} \Omegaop \right)^{-1} h\Omegaop \left[ \fvec(\xvec,t) - \frac{1}{2} \left( \fvec^{\eq}(\xvec,t+h) + \fvec^{\eq}(\xvec,t) \right) \right] \\
%&= \fvec(\xvec,t) - \left( \1 + \frac{h}{2} \Omegaop \right)^{-1} h\Omegaop \left[ \fvec(\xvec,t) - \fvec^{\eq}(\xvec,t) \right] \\
%&= \fvec(\xvec,t) - \Lambdaop \left[ \fvec(\xvec,t) - \fvec^{\eq}(\xvec,t) \right] ,
%\end{split}
%\end{equation}
%%
%}

Equation \eqref{eq:CN-coll} can be rewritten as
%
\begin{equation}
\begin{split}
\mathsf{C} \fvec^\nq &= \left[ \1 - \left( \1 + \frac{h}{2} \Omegaop \right)^{-1} h \Omegaop \right] \, \fvec^\nq \\
&= \left( \1 - \Lambdaop \right) \, \fvec^\nq
\end{split}
\end{equation}
%
with the discrete collision matrix \cite{Dellar2003}
%
\begin{equation}\label{eq:collision-operator}
\Lambdaop = \left( \1 + \frac{h}{2}\Omegaop \right)^{-1}
h \Omegaop .
\end{equation}
%
This corresponds to the familiar replacement of the collision time
$\tau$ by $\tau+h/2$ in the discrete scheme \cite{Henon1987}.
%We will see that this replacement is slightly different in a
%second-order accurate splitting scheme when two half-step collisions
%are applied.
The half-step collision update $\mathsf{C}^\frac12$ in the splitting
scheme \eqref{eq:LB-splitting} is obtained by replacing $h$ by $h/2$
in Eq.~\eqref{eq:collision-operator}
%
\begin{equation}\label{eq:collision-operator2}
\Lambdaop_\frac12 = \left( \1 + \frac{h}{4}\Omegaop \right)^{-1}
\frac{h}{2} \Omegaop .
\end{equation}
%

%\cancel{
%%
%This expression for $\Lambdaop_\frac12$ yields indeed an $O(h^3)$
%approximation of the exact solution (where $\fvec^\eq$ are treated as
%constants again)
%%
%\begin{equation}
%\begin{split}
%\fvec(t+\frac{h}{2}) 
%&= \fvec - ( \1 - e^{-\frac{h}{2}\Omegaop}) \left[ \fvec - \fvec^\eq \right] .\\
%%&= \fvec - \left( \frac{h}{2}\Omegaop - \frac{h^2}{8} \Omegaop^2 \right) \left[ \fvec - \fvec^\eq \right] + O(h^3) \\
%%&= \fvec - \left( \1 + \frac{h}{4} \Omegaop \right)^{-1} \frac{h}{2} \Omegaop \left[ \fvec - \fvec^\eq \right] + O(h^3)\\
%%&= \fvec - \Lambdaop \left[ \fvec - \fvec^\eq \right] + O(h^3) .
%\end{split}
%\end{equation}
%%
%%Note that the formal exact solution rests on the invariance of
%%$\fvec^\eq$ which can this be treated as constants.
%%
%In principle, one could also use the exact expression
%$\1-\exp[-\frac{h}{2}\Omegaop]$ for the discrete collision
%matrix. However, the analysis of Dellar \cite{Dellar2013} revealed
%that the exact solution may lead to spurious decay of shear modes as
%$\tau \lesssim h$. In contrast, the Crank-Nicolson approximation gives
%accurate decay rates even for small collision times. The reason for
%this somewhat counter-intuitive result lies in the compensation of
%errors introduced by the splitting of the Boltzmann equation into
%uncoupled steps on one hand, and the Crank-Nicolson approximation of
%the collision operator on the other.
%
%%Since LBM simulations of soft matter usually operate at low Reynolds
%%number, these considerations are of lesser practical concern in the
%%present context.
%}

%\cancel{
%Using the invariance of $\fvec^\eq$ 
%%
%\begin{equation}
%\begin{split}
%\fvec(t+h) = \mathsf{C} \fvec(t) &= \fvec(t) -  \left( \1 + \frac{h}{2} \Omegaop \right)^{-1} h \Omegaop \left[ \fvec(t) - \fvec^\eq \right] \\
%&= \fvec(t) - \Lambdaop  \left[ \fvec(t) - \fvec^\eq \right] ,
%\end{split}
%\end{equation}
%%
%}

We can convince ourselves that with $\Lambdaop$ and
$\Lambdaop_\frac12$ as in \eqref{eq:collision-operator} and
\eqref{eq:collision-operator2}, two applications of the half-step
collisions $\mathsf{C}^\frac12$ are up to order $O(h^3)$ equivalent to
one full-step collision~$\mathsf{C}$
%
%\cancel{
%\begin{equation}
%\begin{split}
%\mathsf{C}^\frac12\mathsf{C}^\frac12 \fvec
%&= \mathsf{C}^\frac12 \left\{ \fvec^\eq - (\1-\Lambdaop_\frac12) \left[ \fvec - \fvec^\eq \right] \right\} \\
%&= \fvec^\eq + (\1-\Lambdaop_\frac12)^2 \left[ \fvec - \fvec^\eq \right] \\
%&= \fvec^\eq + (\1-\Lambdaop) \left[ \fvec - \fvec^\eq \right] + O(h^3) \\
%&= \mathsf{C} \fvec + O(h^3) .
%\end{split}
%\end{equation}
%}
\begin{equation}
\begin{split}
\mathsf{C}^\frac12\mathsf{C}^\frac12 \fvec^\nq
&= (\1-\Lambdaop_\frac12)^2 \, \fvec^\nq \\
&= (\1-\Lambdaop) \, \fvec^\nq + O(h^3) \\
&= \mathsf{C} \fvec^\nq + O(h^3) .
\end{split}
\end{equation}
%
The order of accuracy follows from $ (\1-\Lambdaop_\frac12)^2 = ( \1 -
\Lambdaop) + O(h^3) $.
%
%\begin{equation}
%\begin{split}
%(1-\Lambdaop_\frac12)^2 %= 1 - \Lambdaop_\frac12 ( 2 - \Lambdaop_\frac12 )
%&= 1 - \left( \1 + \frac{h}{4} \Omegaop \right)^{-2} h \Omegaop \\
%&= 1 - \left( \1 + \frac{h}{2} \Omegaop \right)^{-1} h \Omegaop + O(h^3)
%= ( 1 - \Lambdaop) + O(h^3)
%\end{split}
%\end{equation}
%
Given $\Lambdaop$, we can also exactly define $\1-\Lambdaop_\frac12 =
\sqrt{\1-\Lambdaop}$, where the square-root is to be understood in the
sense of a Cholesky decomposition. For the eigenvalues, we simply have
$1-\lambda_\frac12 = \sqrt{1-\lambda}$.
%
In contrast, the approximation $\mathsf{C}^\frac12 = \frac12(\1 +
\mathsf{C})$ is only $O(h^2)$ accurate since $\Lambdaop_\frac12 =
\Lambdaop/2+O(h^2)$
%
%\cancel{
%\begin{equation}\label{eq:half-col-approx}
%\begin{split}
%\mathsf{C}^\frac12 \fvec 
%&= \fvec - \Lambdaop_\frac12 \left[ \fvec - \fvec^\eq \right] \\
%&= \fvec - \frac{1}{2} \Lambdaop \left[ \fvec - \fvec^\eq \right] + O(h^2)
%= \frac{1}{2} ( \1 + \mathsf{C} ) \fvec + O(h^2) .
%\end{split}
%\end{equation}
%}
\begin{equation}\label{eq:half-col-approx}
\begin{split}
\mathsf{C}^\frac12 \fvec^\nq
&= \left( \1 - \Lambdaop_\frac12 \right) \fvec^\nq \\
&= \left( \1 - \frac{1}{2} \Lambdaop \right) \fvec^\nq + O(h^2)
= \frac{1}{2} ( \1 + \mathsf{C} ) \fvec^\nq + O(h^2) .
\end{split}
\end{equation}
%
When the half-step collisions are applied in every iteration of the
time-marching scheme, the accuracy of $\mathsf{C}^\frac12$ is essential
to maintain the global order $O(h^2)$, cf. section \ref{sec:CSC-LB}.

For the sake of completeness, we mention that in the multi-relaxation
time approach (MRT), the collisions are typically implemented in the
space of moments $m_k$ of the $f_i$:
%
\begin{equation}\label{eq:mrt-coll}
\begin{split}
\mathsf{C} \, m_k^\nq %&= m_k - \lambda_k \left[ m_k - m_k^\eq \right] \\
    &= (1 - \lambda_k) \left[ m_k - m_k^\eq \right] .
\end{split}
\end{equation}
%
Here, $\lambda_k$ is the $k$-th eigenvalue of the
collision matrix which is diagonal in moment space. In the fluctuating
lattice Boltzmann method \cite{Ladd1994a,Adhikari2005,Duenweg2007}, a
noise term is added to the collision step
%
\begin{equation}\label{eq:mrt-fluct}
\mathsf{C} \, m_k^\nq = (1 - \lambda_k) \left[ m_k - m_k^\eq \right] + \sqrt{\rho} \varphi_k r_k ,
\end{equation}
%
where $r_k$ is a Gaussian random number, and the amplitude $\varphi_k
= \sqrt{\mu b_k \lambda_k^2}$ is related to the eigenvalue
$\lambda_k$ and the temperature through $\mu = a^{-3} k_b T /
c_s^2$. $b_k$ is the length of the $k$-th basis vector of the specific
lattice model. Equations \eqref{eq:mrt-coll} and \eqref{eq:mrt-fluct}
are the basis for the implementation of the collision step in
Algorithm~\ref{alg:csc}.  It is worth noting that the fluctuating
collision operator can be split into half-steps by a proper choice of
the eigenvalues and noise strength. A backwards half-step application
$\mathsf{C}^{-\frac12}$ as needed below in equation
\eqref{eq:LB-splitting-combined}, however, is not straightforward to
perform without distorting the correlations of the moments.

%For later reference we introduce the half-step collision operator
%$\mathsf{C}^\frac12$
%%The re-definition of variables for the explicit scheme of He et
%%al. \cite{He} corresponds to applying the collision operator for half
%%a time step. To second order it can be expressed by
%%
%\begin{equation}
%\mathsf{C}^{\frac{1}{2}} = \frac{1}{2} ( \1 + \mathsf{C} ) + O(h^2) ,
%\end{equation}
%%
%%\comment{This is a problem for the overall order!}
%whose application to $\fvec(\xvec,t)$ gives
%%
%\begin{equation}\label{eq:half-collision}
%\begin{split}
%\mathsf{C}^\frac{1}{2} \fvec(\xvec,t) 
%&= \fvec(\xvec,t) - \left( \1 + \frac{h}{2} \Omegaop \right)^{-1} \frac{h}{2} \Omegaop \left[ \fvec(\xvec,t) - \fvec^{\eq}(\xvec,t) \right] \\
%&= \fvec(\xvec,t) - \frac{h}{2} \mathsf{\Lambda} \left[ \fvec(\xvec,t) - \fvec^{\eq}(\xvec,t) \right].
%\end{split}
%\end{equation}
%
%which can be inverted to give
%%
%\begin{equation}
%\bar{\fvec} = \mathsf{C}^{-\frac{1}{2}} \fvec = \fvec + \frac{h}{2}\Omegaop \left[ \fvec - \fvec^\eq \right].
%\end{equation}
%%
%The collision part written in the variables $\bar{\fvec}$ is
%%
%\begin{equation}
%\begin{split}
%\bar{\fvec}(t+h) &= \bar{\fvec}(t) - \left( \1 + \frac{h}{2} \Omegaop
%\right)^{-1} h \Omegaop \left[ \bar{\fvec}(t) - \fvec^\eq \right] \\ 
%&= \bar{\fvec}(t) - h \mathsf{\Lambda} \left[ \bar{\fvec} - \fvec^\eq \right]
%\end{split}
%\end{equation}
%%
%which yields the collision operator in the standard lattice Boltzmann
%algorithm as
%%
%\begin{equation}
%\mathsf{\Lambda} = \left( \1 + \frac{h}{2}\Omegaop \right)^{-1} \Omegaop.
%\end{equation}



\subsection{Forces}
\label{sec:forces}

The force operator $\mathsf{F}$ in Eq.~\eqref{eq:force} can be
discretized using the midpoint rule
%
\begin{subequations}
\begin{equation}\label{eq:force-midpoint}
\begin{split}
\fvec(t+h) &= \mathsf{F} \fvec(t) = \fvec(t) + h \vec{G}(t+\frac{h}{2}) ,
\end{split}
\end{equation}
%
or alternatively with the Crank-Nicolson rule
%
\begin{equation}\label{eq:force-cn}
\begin{split}
\fvec(t+h) &= \mathsf{F} \fvec(t) = \fvec(t) + h \frac{\vec{G}(t+h)+\vec{G}(t)}{2}.
\end{split}
\end{equation}
\end{subequations}
%
The force term $\vec{G}=(G_1,G_2,\dots,G_q)^T$ has the same first
moments as the acceleration term $-\avec \cdot \nabla_\cvec f$
\cite{Luo1998,Guo2002}
%
\begin{align}\label{eq:force-moments}
\sum_i G_i &= 0, & \sum_i \cvec_i G_i &= \Fvec, & \sum_i \cvec_i \cvec_i G_i &= \Fvec \uvec + \uvec \Fvec,
\end{align}
%
where $\Fvec=\rho \avec$. Note that the minus sign for bringing $\avec
\cdot \nabla_\cvec f$ to the right hand side is compensated by another
minus sign that arises in the Hermite expansion of $\nabla_\cvec f$
\cite{Martys1998}.
%
%In a truncated Hermite tensor polynomial expansion it can be expressed
%as \cite{Luo1998,Guo2002}
%%
%\begin{equation}\label{eq:force-hermite}
%G_i = a^{c_i} \left( \frac{\Fvec \cdot \cvec_i}{c_s^2} + \frac{\uvec \cdot (\cvec_i\cvec_i - c_s^2 \1) \cdot \Fvec}{c_s^4} \right) ,
%\end{equation}
%%
%where $a^{c_i}$ are the weights of the lattice and $c_s^2\1=\sum a^{c_i}
%\cvec_i\cvec_i$. 
%
In terms of the first moments, Eqs.~\eqref{eq:force-midpoint} and
\eqref{eq:force-cn} hence read
%
\begin{subequations}\label{eq:LB-force-moments}
\begin{align}
\begin{split}
\uvec(t+h)   &= \uvec(t) + \frac{h}{\rho} \Fvec(t+\frac{h}{2}) ,\\
\Pitens(t+h) &= 
\begin{multlined}[t]
\Pitens(t) + h \big[ \uvec(t+\frac{h}{2}) \Fvec(t+\frac{h}{2})  + \Fvec(t+\frac{h}{2}) \uvec(t+\frac{h}{2}) \big],
\end{multlined}
\end{split}
&& \text{(midpoint)} \\
\intertext{and}
\begin{split}
\uvec(t+h) &= \uvec(t) + \frac{h}{2\rho} \left( \Fvec(t+h) + \Fvec(t)
\right),\\
\Pitens(t+h) &= 
\begin{multlined}[t][.5\textwidth]
\Pitens(t) + \frac{h}{2} \big[ \uvec(t) \Fvec(t) +
\Fvec(t) \uvec(t) \\ + \uvec(t+h) \Fvec(t+h) + \Fvec(t+h) \uvec(t+h)
\big] .
\end{multlined}
\end{split}
&& \text{(Crank-Nicolson)}
\end{align}
\end{subequations}
%
If the force is constant during a time step $h$, Crank-Nicolson and
midpoint are both equivalent to forward Euler discretization, where
the midpoint velocity $\uvec(t+h/2)$ has to be used in the stress
update.
%

%Applying the splitting scheme for $N$ timesteps reads
%\begin{equation}
%\begin{split}
%f_i(\{\xvec\},t+Nh) &= \left[ \mathsf{F}^{\frac{1}{2}} \mathsf{C}^{\frac{1}{2}} \mathsf{S} \mathsf{C}^{\frac{1}{2}} \mathsf{F}^{\frac{1}{2}} \right] ^N f_i(\{\xvec\},t) \\
%&= \mathsf{F}^\frac{1}{2} \left[ \mathsf{C}^\frac{1}{2} \mathsf{S} \mathsf{C}^\frac{1}{2} \mathsf{F} \right]^N \mathsf{F}^{-\frac{1}{2}} f_i(\{\xvec\},t) \\
%&= \mathsf{F}^\frac{1}{2} \mathsf{C}^\frac{1}{2} \left[ \mathsf{S} \mathsf{C}^\frac{1}{2} \mathsf{F} \mathsf{C}^\frac{1}{2} \right]^N \mathsf{C}^{-\frac{1}{2}} \mathsf{F}^{-\frac{1}{2}} f_i(\{\xvec\},t)
%\end{split}
%\end{equation}

It is worthwhile to note that up to $O(\Fvec^2)$, the force operator
$\mathsf{F}$ only acts on the equilibrium distribution and leaves the
non-equilibrium part unchanged
%
\begin{equation}
\begin{split}
%f_i^\eq(\rho,\uvec(t+h)) &= f_i^\eq(\rho,\uvec(t)) + h G_i(t+\frac{h}{2}) + O(\Fvec^2)\\
\fvec(t+h) - \fvec^\eq(\rho,\uvec(t+h)) &= \fvec(t) + h \vec{G}(t+\frac{h}{2}) - \fvec^\eq(\rho,\uvec(t+h)) \\
&= \fvec(t) - \fvec^\eq(\rho,\uvec(t)) + O(\Fvec^2) ,
\end{split}
\end{equation}
%
which can easily be checked using the truncated Hermite expansion of
the force term $G_i$ and the equilibrium distribution $f_i^\eq$
\cite{Luo1998,Martys1998,Guo2002}.
%
%\begin{equation}
%\begin{split}
%&f_i(t+h) = f_i^\eq(\rho,\uvec(t)) + \left[ f_i(t) - f_i^\eq(\rho,\uvec(t)) \right] + h G_i(t+h/2)\\
%&= 
%\begin{multlined}[t]
%a^{c_i} \rho \left( 1 + \frac{\uvec(t)\cdot\cvec_i}{c_s^2} 
%+ \frac{\uvec(t) \cdot (\cvec_i\cvec_i - c_s^2 \1) \cdot \uvec(t)}{2c_s^4} \right) + \left[ f_i(t) - f_i^\eq((\rho,\uvec(t)) \right] \\
%+ h a^{c_i} \left( \frac{\Fvec(t+h/2) \cdot \cvec_i}{c_s^2} 
%+ \frac{\uvec(t+h/2) \cdot (\cvec_i\cvec_i - c_s^2 \1) \cdot \Fvec(t+h/2)}{c_s^4} \right) \\
%\end{multlined} \\
%&= 
%\begin{multlined}[t]
%a^{c_i} \rho \left( 1 + \frac{\uvec(t+h)\cdot\cvec_i}{c_s^2} 
%+ \frac{\uvec(t+h) \cdot (\cvec_i\cvec_i - c_s^2 \1) \cdot \uvec(t+h)}{2 c_s^4} \right) + O(\Fvec^2)\\
%+ \left[ f_i(t) - f_i^\eq((\rho,\uvec(t)) \right] 
%\end{multlined} \\
%&= f_i^\eq(\rho,\uvec(t+h)) + \left[ f_i(t) - f_i^\eq((\rho,\uvec(t)) \right] + O(\Fvec^2)
%\end{split}
%\end{equation}
%
As a consequence, the force and collision operators approximately
commute
%
\begin{equation}
\begin{split}
\mathsf{C}\mathsf{F} \fvec(t) &= \mathsf{F}\fvec(t) - \Lambdaop \mathsf{F} \left[ \fvec(t) - \fvec^\eq(\rho,\uvec(t)) \right] \\
&= \fvec(t) + h \vec{G} - \Lambdaop \left[ \fvec(t) - \fvec^\eq(\rho,\uvec(t)) \right] + O(\Fvec^2) \\
&= \mathsf{C} \fvec(t) + h \vec{G} + O(\Fvec^2) \\
&= \mathsf{F}\mathsf{C} \fvec(t) + O(\Fvec^2) .
\end{split}
\end{equation}
%
For small enough forces the order of separate collisions and forcing
thus does not matter. Note that throughout this section we have
discussed a full time step. The action of the half-step force operator
$\mathsf{F}^\frac12$ is obtained by simply replacing $h$ by $h/2$. We
will see below that the splitting scheme leads to the familiar
transformation of the force term required for second-order accuracy
\cite{Guo2002,He1998a,Luo1998}.



\section{Collide-stream-collide lattice Boltzmann}
\label{sec:CSC-LB}

Strang or Trotter-Suzuki splitting as in Eq.~\eqref{eq:LB-splitting}
leads to a sequence of operations which involve the half-step collision
operator $\mathsf{C}^\frac12$. If it is the first and last operation,
then upon concatenation of time-steps the adjacent collision
half-steps can be combined, and the standard LB algorithm is recovered
\cite{Dellar2013}
%
\begin{equation}\label{eq:LB-splitting-combined}
\left[ \mathsf{C}^\frac12 \mathsf{S} \mathsf{C}^\frac12 \right] ^n = \
\mathsf{C}^\frac12 \left[ \mathsf{S} \mathsf{C} \right]^n
\mathsf{C}^{-\frac12} ,
\end{equation}
%
where the forces are included in the collision operator. At the
beginning and end of the simulation the half-step collisions
$\mathsf{C}^{-\frac12}$ and $\mathsf{C}^{\frac12}$ have to be applied,
respectively.
%
In this work, we group the operations such that the force operator
$\mathsf{F}$ does not interrupt the collision and streaming steps,
%
%the action of the forces is applied in separate
%steps. This leaves us with two options for the splitting (as mentioned
%above the streaming operator must not be split).
%\begin{align}
%\left[  \mathsf{C}^{\frac{1}{2}} \mathsf{F}^{\frac{1}{2}} \mathsf{S} \mathsf{F}^{\frac{1}{2}} \mathsf{C}^{\frac{1}{2}} \right]^N
%&= \mathsf{C}^\frac{1}{2} \mathsf{F}^\frac{1}{2} \left[ \mathsf{S} \mathsf{F}^\frac{1}{2} \mathsf{C} \mathsf{F}^\frac{1}{2} \right]^N \mathsf{F}^{-\frac{1}{2}} \mathsf{C}^{-\frac{1}{2}} f_i(\{\xvec\},t) \\
%\left[  \mathsf{F}^{\frac{1}{2}} \mathsf{C}^{\frac{1}{2}} \mathsf{S} \mathsf{C}^{\frac{1}{2}} \mathsf{F}^{\frac{1}{2}} \right]^N
%&= \mathsf{F}^\frac{1}{2} \mathsf{C}^\frac{1}{2} \left[ \mathsf{S} \mathsf{C}^\frac{1}{2} \mathsf{F} \mathsf{C}^\frac{1}{2} \right]^N \mathsf{C}^{-\frac{1}{2}} \mathsf{F}^{-\frac{1}{2}} f_i(\{\xvec\},t) \\
%&= \mathsf{F}^\frac{1}{2} \left[ \mathsf{C}^\frac{1}{2} \mathsf{S} \mathsf{C}^\frac{1}{2} \mathsf{F} \right]^N \mathsf{F}^{-\frac{1}{2}} f_i(\{\xvec\},t) 
%\end{align}
%In the last equation, the operators are grouped such that the force
%operator $\mathsf{F}$ does not interrupt the collision and streaming
%process. 
%
which facilitates the implementation of hybrid coupling algorithms with
sub-splitting of $\mathsf{F}$ (cf. section~\ref{sec:multi-timestep}).

Before developing the details of the coupling, we analyze the
splitting $\mathsf{C}^\frac12 \mathsf{S} \mathsf{C}^\frac12$ for the
LB part. The latter suggests a collide-stream-collide lattice
Boltzmann update with two half-step collisions before and after the
streaming. Let $f_i^*$ denote the populations after the first
collision step, and $\bar f_i$ the populations after the streaming
step. Then the three separate steps are
%
%What needs to be clarified though is the notion of a collision
%half-step $\mathsf{C}^\frac{1}{2}$. \comment{Accuracy?}

%Recall that $\mathsf{C}^\frac12=\frac{1}{2}(\1+\mathsf{C})$ and
%\begin{equation}
%\begin{split}
%\mathsf{C}^\frac12 \fvec &= \fvec - \frac{1}{2} \left(\1 + \frac{h}{2}\Omegaop \right)^{-1} h \Omegaop \left[ \fvec - \fvec^\eq \right]\\
%&= \fvec -\frac{h}{2} \mathsf{\Lambda} \left[ \fvec - \fvec^\eq \right].
%\end{split}
%\end{equation}

%As can be seen from Eq.~\eqref{eq:half-collision}, a collision
%half-step relaxes the non-equilibrium part towards the equilibrium at
%half the relaxation rate. The collide-stream-collide lattice Boltzmann
%update is thus
%%
%\begin{equation}
%\begin{split}
%f_i(\xvec+h\cvec_i,t+h) &= f'_i(\xvec+h\cvec_i,t+h) - \frac{h}{2} \sum_j \mathsf{\Lambda}_{ij} \left[ f'_j(\xvec+h\cvec_i,t+h) - f_j^\eq(\xvec+h\cvec_i,t+h) \right] \\
%&=f_i^*(\xvec,t) - \frac{h}{2} \sum_j \mathsf{\Lambda}_{ij} \left[ f_j^*(\xvec,t) - f_j^\eq(\xvec+h\cvec_i,t+h) \right] \\
%&= f_i(\xvec,t) - \frac{h}{2} \sum_j \mathsf{\Lambda}_{ij} \left[ f_j(\xvec,t) - f_j^\eq(\xvec,t) \right] \\
%&\qquad - \frac{h}{2} \sum_j \mathsf{\Lambda}_{ij} \left\{ f_j(\xvec,t) - \frac{h}{2} \sum_k \mathsf{\Lambda}_{jk} \left[ f_k(\xvec,t) - f_k^\eq(\xvec,t) \right] - f_j^\eq(\xvec+h\cvec_i,t+h) \right\} \\
%&= f_i(\xvec,t) - \frac{h}{2} \sum_j \mathsf{\Lambda}_{ij} \left[ f_j(\xvec,t) - f_j^\eq(\xvec,t) \right] - \frac{h}{2} \sum_j \mathsf{\Lambda}_{ij} \left[ f_j(\xvec,t) - f_j^\eq(\xvec+h\cvec_i,t+h) \right] 
%+ O(h^2) \\
%&= f_i(\xvec,t) - h \sum_j \mathsf{\Lambda}_{ij} \left[ f_j(\xvec,t) - \frac{f_j^\eq(\xvec+h\cvec_i,t+h) + f_j^\eq(\xvec,t)}{2} \right] + O(h^2) \\
%&= f_i(\xvec,t) - h \sum_j \mathsf{\Lambda}_{ij} \left[ f_j(\xvec,t) - f_j^\eq(\xvec,t) \right] + O(h^2) ,
%\end{split}
%\end{equation}
%where $f_i^*$ denotes the populations after the first collision step,
%and $f_i'$ the populations after the streaming step. Up to $O(h^2)$
%errors the split scheme
%$\mathsf{C}^\frac12\mathsf{S}\mathsf{C}^\frac12$ is equivalent to the
%standard LB collide-stream $\mathsf{S}\mathsf{C}$, hence the
%derivation of macroscopic equations by the Chapman-Enskog expansion
%\cite{ChapmanCowling} remains the same. The relation of the collision
%rates to the viscosity of the fluid is thus known
%\cite{McNamara1993,Ladd2001,Dellar2001}. The application of a
%collision half-step is thus easily achieved by multiplying the
%relaxation matrix with the factor $\frac{1}{2}$ without affecting
%overall second-order accuracy. \comment{This is probably not true!}
%
%The collide-stream-collide lattice Boltzmann update with half-step
%collision matrix $\Lambdaop_\frac12$ is given by the steps
%
\begin{subequations}\label{eq:csc}
\begin{align}
&\fvec(\xvec,t) = 
\begin{multlined}[t]
\bar \fvec(\xvec,t) - \Lambdaop_{\frac12} \left[ \bar \fvec(\xvec,t) - \fvec^\eq(\xvec,t) \right]
\end{multlined} \label{eq:trafoHe} \\
&\bar f_i(\xvec+h\cvec_i,t+h) = f_i^*(\xvec,t) \label{eq:poststreaming} \\
&\fvec^*(\xvec,t) =
\begin{multlined}[t]
\fvec(\xvec,t) - \Lambdaop_\frac12 \left[ \fvec(\xvec,t) - \fvec^\eq(\xvec,t) \right]
\end{multlined}.
%\intertext{or combined}
%&\fvec(\xvec+h\cvec,t+h) =
%\begin{multlined}[t]
%\fvec(\xvec,t) - \Lambdaop_\frac12 (2 - \Lambdaop_\frac12) \left[ \fvec(\xvec,t) - \fvec^\eq(\xvec,t) \right] \\
%+ \Lambdaop_\frac12 \left[ \fvec^\eq(\xvec+h\cvec,t+h) - \fvec^\eq(\xvec,t) \right] .
%\end{multlined}
\end{align}
\end{subequations}
%
%\begin{equation}
%\begin{split}
%\fvec(\xvec+h\cvec,t+h) &= 
%\begin{multlined}[t]
%\fvec'(\xvec+h\cvec,t+h) \\- \Lambdaop_{\frac12} \left[ \fvec'(\xvec+h\cvec,t+h) - \fvec^\eq(\xvec+h\cvec,t+h) \right]
%\end{multlined}\\
%&= \fvec^*(\xvec,t) - \Lambdaop_\frac12 \left[ \fvec^*(\xvec,t) - \fvec^\eq(\xvec+h\cvec,t+h) \right] \\
%&=
%\begin{aligned}[t]%multlined}[t][.67\textwidth]
%\fvec(\xvec,t) &- \Lambdaop_\frac12 \left[ \fvec(\xvec,t) - \fvec^\eq(\xvec,t) \right] \\
%&
%\begin{multlined}[t][.59\textwidth]
%- \Lambdaop_\frac12 \left[ \fvec(\xvec,t)  - \Lambdaop_\frac12 \left( \fvec(\xvec,t) - \fvec^\eq(\xvec,t) \right) \right. \\
%- \fvec^\eq(\xvec+h\cvec,t+h) \Bigr]
%\end{multlined}
%\end{aligned}\\
%%&= \left( 1 - \Lambda_\frac12 \right) \fvec'(\xvec+h\cvec,t+h) + \Lambda^\frac12 \fvec^\eq(\xvec+h\cvec,t+h) \\
%%&= \left( 1 - \Lambda_\frac12 \right) \fvec^*(\xvec,t) + \Lambdaop_\frac12 \fvec^\eq(\xvec+h\cvec,t+h) \\
%%&= \left( 1 - \Lambdaop_\frac12 \right) \left[ \fvec^\eq(\xvec,t) + \left( 1 - \Lambdaop_\frac12 \right) \left( \fvec(\xvec,t) - \fvec^\eq(\xvec,t) \right) \right] + \Lambdaop_\frac12 \fvec^\eq(\xvec+h\cvec,t+h) \\
%%&= \fvec^\eq(\xvec,t) + \left( 1 - \Lambdaop_\frac12 \right)^2 \left[ \fvec(\xvec,t) - \fvec^\eq(\xvec,t) \right] + \Lambdaop_\frac12 \left[ \fvec^\eq(\xvec+h\cvec,t+h) - \fvec^\eq(\xvec,t) \right] \\
%&=
%\begin{aligned}[t]%{multlined}[t]
%\fvec(\xvec,t) &- 2 \Lambdaop_\frac12 \left[ \fvec(\xvec,t) - \frac{\fvec^\eq(\xvec+h\cvec,t+h) + \fvec^\eq(\xvec,t)}{2} \right] \\
%&+ \Lambdaop_\frac12 \Lambdaop_\frac12 \left[ \fvec(\xvec,t) - \fvec^\eq(\xvec,t) \right] 
%\end{aligned}\\%{multlined}\\
%&=
%\begin{aligned}[t]
%\fvec(\xvec,t) &- \Lambdaop_\frac12 (2 - \Lambdaop_\frac12) \left[ \fvec(\xvec,t) - \fvec^\eq(\xvec,t) \right] \\
%&+ \Lambdaop_\frac12 \left[ \fvec^\eq(\xvec+h\cvec,t+h) - \fvec^\eq(\xvec,t) \right] , 
%\end{aligned}\\
%%&= \fvec(\xvec,t) - \left( 1 - (1 - \Lambdaop_\frac12)^2 \right) \left[ \fvec(\xvec,t) - \fvec^\eq(\xvec,t) \right] + \Lambdaop_\frac12 \left[ \fvec^\eq(\xvec+h\cvec,t+h) - \fvec^\eq(\xvec,t) \right] 
%\end{split}
%\end{equation}
%
The difference to a direct Crank-Nicolson discretization of the
Boltzmann equation is that two explicit collisions are applied, and
not just one collision using a midpoint approximation for the
non-equilibrium distribution~\cite{Xu2003,Ubertini2010}.
%
%The first line of Eq.~\eqref{eq:csc} can be inverted to give
%%
%\begin{equation}\label{eq:trafoHe}
%\begin{split}
%\fvec(\xvec,t) &= 
%%\fvec'(\xvec,t) - \Lambdaop_\frac12 \left[ \fvec'(\xvec,t) - \fvec^\eq(\xvec,t) \right] = 
%\fvec^\eq(\xvec,t) + (1-\Lambdaop_\frac12) \left[ \fvec'(\xvec,t) - \fvec^\eq(\xvec,t) \right] \\
%\fvec'(\xvec,t) &= 
%%\fvec(\xvec,t) + ( 1 - \Lambdaop_\frac12 )^{-1} \Lambdaop_\frac12 \left[ \fvec(\xvec,t) - \fvec^\eq(\xvec,t) \right] = 
%\fvec^\eq(\xvec,t) + (1-\Lambdaop_\frac12)^{-1} \left[ \fvec(\xvec,t) - \fvec^\eq(\xvec,t) \right] . \\
%\end{split}
%\end{equation}
%
To see the relation to the standard LB scheme, we note that equation
\eqref{eq:poststreaming} for the post-streaming distribution can be
rewritten as
%
\begin{equation}
\begin{split}
\bar f_i(\xvec+h\cvec_i,t+h) &= f_i^*(\xvec,t) \\ 
%&= \fvec(\xvec,t) - \Lambdaop_\frac12 \left[ \fvec(\xvec,t) - \fvec^\eq(\xvec,t) \right] \\
&= f_i^\eq(\xvec,t) + \sum_j (\1-\Lambdaop_\frac12)_{ij} \left[ f_j(\xvec,t) -
f_j^\eq(\xvec,t) \right] \\ 
&= f_i^\eq(\xvec,t) + \sum_j (\1-\Lambdaop_\frac12)^2_{ij} \left[ \bar f_j(\xvec,t) - f_j^\eq(\xvec,t)
\right] , \\
%&= ( 1 - \Lambdaop_\frac12 ) \fvec(\xvec,t) + \Lambdaop_\frac12 \fvec^\eq(\xvec,t) \\
%&= ( 1 - \Lambdaop_\frac12 ) \left[ ( 1 - \Lambdaop_\frac12 ) \fvec'(\xvec,t) + \Lambdaop_\frac12 \fvec^\eq(\xvec,t) \right] + \Lambdaop_\frac12 \fvec^\eq(\xvec,t) \\
%&= \left( 1 - (2 - \Lambdaop_\frac12) \Lambdaop_\frac12 \right) \fvec'(\xvec,t) + ( 2 - \Lambdaop_\frac12 ) \Lambdaop_\frac12 \fvec^\eq(\xvec,t) \\
%&= \fvec'(\xvec,t) - (2 - \Lambdaop_\frac12) \Lambdaop_\frac12 \left[ \fvec'(\xvec,t) - \fvec^\eq(\xvec,t) \right] \\
%\fvec^*(\xvec,t) &= \fvec(\xvec,t) - \Lambdaop_\frac12 \left[ \fvec(\xvec,t) - \fvec^\eq(\xvec,t) \right] \\
%\fvec(\xvec,t) &= \fvec^*(\xvec, t) + ( 1 - \Lambdaop_\frac12 )^{-1} \Lambdaop_\frac12 \left[ \fvec^*(\xvec,t) - \fvec^\eq(\xvec,t) \right]
\end{split}
\end{equation}
%
which is the standard LB scheme with collision matrix
$(1-\Lambdaop_\frac12)^2$. A Chapman-Enskog expansion
\cite{ChapmanCowling} reveals that the moments of $\bar f_i$ satisfy
the Navier-Stokes equation, which allows us to obtain explicit
expressions for the transport coefficients.
%
%On the other hand, it follows from the second-order accuracy of the
%splitting scheme and Eq.~\eqref{eq:trafoHe} that
%%
%\begin{equation}
%\tens{\sigma} = \sum_i (f_i - f^\eq_i) \cvec_i \cvec_i = (1-\Lambdaop_\frac12) \sum_i (f'_i - f^\eq_i) \cvec_i \cvec_i .
%\end{equation}
%%
%The difference between these two expressions is of order $O(h^2)$
%which maintains overall second-order accuracy. 
%
%Eq.~\eqref{eq:csc} can be rearranged
%%
%\begin{equation}
%\begin{split}
%\fvec(\xvec+h\cvec,t+h) - \Lambdaop_\frac12 \fvec^\eq(\xvec+h\cvec,t+h) &= \fvec(\xvec,t) - \Lambdaop_\frac12 \fvec^\eq(\xvec,t) - \Lambdaop_\frac12 (2-\Lambdaop_\frac12) \left[ \fvec(\xvec,t) - \fvec^\eq(\xvec,t) \right] \\
%\tilde\fvec(\xvec+h\cvec,t+h) &= \tilde\fvec(\xvec,t) - \Lambdaop_\frac12 (2-\Lambdaop_\frac12) \left[ \tilde\fvec(\xvec,t) - (1-\Lambdaop_\frac12) \fvec^\eq(\xvec,t) \right] \\
%&= \tilde\fvec(\xvec,t) - \left( 1 - ( 1 - \Lambdaop_\frac12 )^2 \right) \left[ \tilde\fvec(\xvec,t) - (1-\Lambdaop_\frac12) \fvec^\eq(\xvec,t) \right] 
%\end{split}
%\end{equation}
%%
%where we have introduced the transformation
%%
%\begin{equation}
%\tilde{\fvec} = \fvec - \Lambdaop_\frac12 \fvec^\eq .
%\end{equation}
%%
%The standard LB scheme is recovered if $( 1 - \Lambdaop_\frac12 ) =
%\sqrt{1-\Lambda}$, because
%%
%\begin{equation}
%\mathsf{C}^\frac12 \mathsf{C}^\frac12 \fvec = \fvec - \Lambdaop_\frac12 \left( 2 - \Lambdaop_\frac12 \right) \left[ \fvec - \fvec^\eq \right] = \fvec - \left( 1 - ( 1 - \Lambdaop_\frac12 )^2 \right) \left[ \fvec - \fvec^\eq \right] = \fvec - \Lambdaop \left[ \fvec - \fvec^\eq \right] = \mathsf{C} \fvec .
%\end{equation}
%%
%The $\tilde\fvec$ satisfy Navier-Stokes per Chapman-Enskog.
%
%The $\fvec$ satisfy the Boltzmann equation.
%
%Stress/viscosity from the two should be different.
%
%\begin{equation}
%\tilde m_k = m - \lambda_k m_k^\eq
%\end{equation}
%
%conserved moments $\gamma_k = 1$
%
%$1-\lambda_k=\gamma_k$
%
%\begin{equation}
%\tilde{\overline{\tens{\Pi}}} = \overline{\tens{\Pi}} - \lambda_s \overline{\tens{\Pi}}^\eq
%\end{equation}
%
%The collide-stream-collide scheme is yet slightly different, as we
%will see in the following asymptotic expansion. A Taylor expansion of
%the last line in Eq.~\eqref{eq:csc} leads to
%%
%\begin{equation}\label{eq:csc-taylor}
%\begin{split}
%\left[ h \left( \frac{\d}{\d t} + \cvec_i  \frac{\d}{\d \xvec} \right)
%+\frac{h^2}{2} \left( \frac{\d}{\d t} + \cvec_i  \frac{\d}{\d \xvec} \right) ^ 2 \right]
%\left( f_i - \Lambdaop_\frac12 f_j^\eq \right) 
%&= - \Lambdaop_\frac12 (2 - \Lambdaop_\frac12) \left[ f_k - f_k^\eq \right] ,
%\end{split}
%\end{equation}
%%
%where we have grouped all derivative terms on the left hand side. Note
%that the right hand side gives the collisional contribution for a full
%time step if $( 1 - \Lambdaop_\frac12 ) = \sqrt{1-\Lambda}$, because
%%
%\begin{equation}
%\mathsf{C}^\frac12 \mathsf{C}^\frac12 \fvec = \fvec - \Lambdaop_\frac12 \left( 2 - \Lambdaop_\frac12 \right) \left[ \fvec - \fvec^\eq \right] = \fvec - \left( 1 - ( 1 - \Lambdaop_\frac12 )^2 \right) \left[ \fvec - \fvec^\eq \right] = \fvec - \Lambdaop \left[ \fvec - \fvec^\eq \right] = \mathsf{C} \fvec .
%\end{equation}
%The square root $\sqrt{1-\Lambda}$ is understood in the sense of a
%Cholesky decomposition.
%
%A Chapman-Enskog expansion of equation \eqref{eq:csc-taylor} leads to
%the following hierarchy of equations:
%%
%\begin{equation}
%\begin{split}
%0 &= - \Lambdaop_\frac12 (2 - \Lambdaop_\frac12) \left[ f_k^{(0)} -
%f_k^\eq \right] \\
%h \left( \frac{\d}{\d t_1} + e_{\gamma i} \frac{\d}{\d x_\gamma} \right)
%\left( f_i^{(0)} - \Lambdaop_\frac12 f_j^\eq \right) &= -
%\Lambdaop_\frac12 (2 - \Lambdaop_\frac12) \ f_k^{(1)} \\
%%\left[ h \frac{\d}{\d t_2}  +\frac{h^2}{2} \left( \frac{\d}{\d t} + \cvec_i  \frac{\d}{\d \xvec} \right) ^ 2 \right] \left( f_i^{(0)} - \Lambdaop_\frac12 f_i^\eq \right)
%%+ h \left( \frac{\d}{\d t_1} + \cvec_i \frac{\d}{\d \xvec} \right) \ f_i^{(1)}
%%&= - \Lambdaop_\frac12 (2 - \Lambdaop_\frac12) \ f_i^{(2)} \\
%h \frac{\d}{\d t_2} \left( f_i^{(0)} - \Lambdaop_\frac12 f_j^\eq \right) 
%+\frac{h}{2} \left( \frac{\d}{\d t} + e_{\gamma i}  \frac{\d}{\d x_\gamma} \right)
%\left( 2  - \Lambdaop_\frac12 (2 - \Lambdaop_\frac12) \right) f_k^{(1)}
%&= - \Lambdaop_\frac12 (2 - \Lambdaop_\frac12) \ f_k^{(2)} \\
%\end{split}
%\end{equation}
%%
%The difference to the Chapman-Enskog expansion of the standard LB
%scheme are the terms $f_i^{(0)}- \Lambdaop_\frac12 f_j^\eq$ where
%normally only $f_i^{(0)}$ would appear.
%
%Can these terms be absorbed in some redifinition $\tilde{f}_i$?
%
%The moment equations on the $t_1$-scale are
%%
%\begin{equation}
%\begin{split}
%\d_{t_1} \rho + \d_\gamma j_\gamma &= 0 \\
%\d_{t_1} j_\alpha + \d_\beta ( 1 - \lambda_e ) \Pitens^{(0)}_{\alpha\beta} &= 0 \\
%\d_{t_1} ( 1 - \lambda_e) \Pitens^{(0)}_{\alpha\beta} + \d_\gamma ( 1 - \lambda_o ) \tens{\Phi}^{(0)}_{\alpha\beta\gamma} &= - \Lambdaop_\frac12 ( 2 - \Lambdaop_\frac12 ) \Pitens^{(1)}_{\alpha\beta}
%\end{split}
%\end{equation}
%
%\begin{equation}
%\begin{split}
%&\phantom{=} - ( 1 - \lambda_e ) c_s^2 \left[ \delta_{\alpha\beta} \d_\gamma j_\gamma + ( 1 - \lambda_e ) \left( u_\alpha \d_\beta \rho + u_\beta \d_\alpha \rho \right) \right] + ( 1 - \lambda_o) c_s^2 \left[ \delta_{\alpha\beta} \d_\gamma j_\gamma + \d_\alpha j_\beta + \d_\beta j_\alpha \right]\\
%&= c_s^2 (\lambda_e - \lambda_o) \delta_{\alpha\beta} \d_\gamma j_\gamma 
%   + c_s^2 \left[ (1-\lambda_o) \d_\alpha j_\beta - (1-\lambda_e)^2 u_\beta \d_\alpha \rho \right]
%   + c_s^2 \left[ (1-\lambda_o) \d_\beta j_\alpha - (1-\lambda_e)^2 u_\alpha \d_\beta \rho \right] \\
%&= c_s^2 (\lambda_e - \lambda_o) \delta_{\alpha\beta} \d_\gamma j_\gamma 
%   + c_s^2 \left[ \rho \d_\alpha u_\beta - \lambda_o \d_\alpha j_\beta + \lambda_e (2 - \lambda_e) u_\beta \d_\alpha \rho \right]
%   + c_s^2 \left[ \rho \d_\beta u_\alpha - \lambda_o \d_\beta j_\alpha + \lambda_e (2 - \lambda_e) u_\alpha \d_\beta \rho \right]
%\end{split}
%\end{equation}

The viscosity $\nu$ is given
by~\cite{McNamara1993,Dellar2001,Ladd2001,Duenweg2008}
%
\begin{equation}\label{eq:viscosity}
\nu = \frac{h c_s^2}{2} \frac{1+\gamma_\frac12^2}{1-\gamma_\frac12^2} 
%= \frac{h c_s^2}{2} \frac{2 - 2 \lambda_\frac12 + \lambda_\frac12^2}{2\lambda_\frac12 - \lambda_\frac12^2}  
= {h c_s^2} \left( \frac{1}{2\lambda_\frac12-\lambda_\frac12^2} - \frac{1}{2} \right)
= {h c_s^2} \left( \frac{1}{2\lambda_\frac12} - \frac{1}{4} \right) +O(h^2) ,
%= {h c_s^2} \left( \frac{1}{\lambda} - \frac{1}{2} \right)
\end{equation}
%
where $\gamma_\frac12 = 1 - \lambda_\frac12$ is the eigenvalue of $( 1
- \Lambdaop_\frac12)$ corresponding to the shear modes. The relation
to the continuous Boltzmann equation with BGK collisions
\cite{Bhatnagar1954} is achieved by shifting the BGK relaxation time
$\tau$ according to the time discretization~\cite{Henon1987,Xu2003,
  Ubertini2010}, cf. Eq.~\eqref{eq:collision-operator2},
%
\begin{equation}\label{eq:BGKlambda}
\lambda_\frac12 = \frac{h}{2} \frac{1}{\tau + \frac{h}{4}} ,
\end{equation}
%
such that
%
\begin{equation}
\nu = {h c_s^2} \left(\frac{\tau+\frac{h}{4}}{h} -  \frac{1}{4} \right) = \tau c_s^2 ,
%\frac{h c_s^2}{2} \frac{1 + (1-\lambda_\frac12)^2}{1 - (1-\lambda_\frac12)^2
%= \frac{h c_s^2}{2} \frac{\left(\tau+\frac{h}{4}\right)^2 + \left(\tau-\frac{h}{4}\right)^2}{\left(\tau+\frac{h}{4}\right)^2 - \left(\tau-\frac{h}{4}\right)^2} =  \frac{h c_s^2}{2} \frac{2\tau^2 + \frac{h^2}{8}}{h \tau} = \tau c_s^2 
%+ O\left(h^2\right)
\end{equation}
%
in accordance with the continuum result \cite{ChapmanCowling}.

The viscous stress $\tens{\sigma}$ is given
by~\cite{Duenweg2008,Dellar2001}
%
\begin{equation}\label{eq:viscstress}
\begin{split}
\tens{\sigma} &= \frac{1}{2} \left( \1 + (\1-\Lambdaop_\frac12)^2 \right) \sum_i (\bar f_i - f^\eq_i) \cvec_i \cvec_i \\
&= \frac{(\1-\Lambdaop_\frac12)^{-1} + (\1-\Lambdaop_\frac12)}{2} \sum_i (f_i - f^\eq_i) \cvec_i \cvec_i \\
&= \sum_i (f_i - f^\eq_i) \cvec_i \cvec_i + O(h^2) .
\end{split}
\end{equation}
%
The second line corresponds to a midpoint scheme that averages the
stress before the last half-step collision and after the next
half-step collision. The last line shows that in the
collide-stream-collide scheme, the viscous stress is equal to the
local non-equilibrium stress up to $O(h^2)$ and can directly be
obtained from the distributions $f_i$. No additional transformation or
look-ahead collision as in the standard scheme has to be performed.

%Comparing the half-step collisions to the full-step collisions of the
%standard LB scheme, we find \comment{Check and compare with above and clean up}
%%
%\begin{equation}
%\mathsf{C}^\frac12 \mathsf{C}^\frac12 \fvec = \fvec^\eq + \left( 1 - \Lambdaop_\frac12 \right)^2 \left[ \fvec - \fvec^\eq \right] = \fvec^\eq +(1 - \Lambdaop) \left[ \fvec - \fvec^\eq \right] = \mathsf{C} \fvec ,
%\end{equation}
%%
%from which it follows $1-\Lambdaop_\frac12=\sqrt{1-\Lambdaop}$. The
%square-root is to be understood in the sense of a Cholesky
%decomposition. For the eigenvalues, we simply have $\gamma_\frac12 = 1
%- \lambda_\frac12 = \sqrt{1 - \lambda} = \sqrt{\gamma}$.
%
%\begin{equation}
%\begin{split}
%1 - \Lambdaop_\frac12 &= 1 - \left( \mathsf{I} + \frac{h}{4} \mathsf{\Omega} \right)^{-1} \frac{h}{2} \mathsf{\Omega} = \left( \mathsf{I} + \frac{h}{4} \mathsf{\Omega} \right)^{-1} \left( \mathsf{I} - \frac{h}{4} \mathsf{\Omega} \right) 
%= \left( \mathsf{I} - \frac{h}{4} \mathsf{\Omega} \right)^2 + O(h^2) \\
%( 1 - \Lambdaop_\frac12 ) ^ 2 &= \left( \mathsf{I} + \frac{h}{4} \mathsf{\Omega} \right)^{-1} \left( \mathsf{I} - \frac{h}{4} \mathsf{\Omega} \right) \left( \mathsf{I} + \frac{h}{4} \mathsf{\Omega} \right)^{-1} \left( \mathsf{I} - \frac{h}{4} \mathsf{\Omega} \right) \\
%&= \left( \mathsf{I} - \left( \mathsf{I} + \frac{h}{4} \mathsf{\Omega} \right)^{-1} \frac{h}{2} \mathsf{\Omega} \right) \left( \mathsf{I} + \frac{h}{4} \mathsf{\Omega} \right)^{-1} \left( \mathsf{I} - \frac{h}{4} \mathsf{\Omega} \right) \\
%&= \left( \mathsf{I} + \frac{h}{4} \mathsf{\Omega} \right)^{-1} \left( \mathsf{I} - \frac{h}{2}\mathsf{\Omega} \left(  \mathsf{I} + \frac{h}{4} \mathsf{\Omega} \right)^{-1} \right) \left( \mathsf{I} - \frac{h}{4} \mathsf{\Omega} \right) \\
%&= \left( \mathsf{I} + \frac{h}{4} \mathsf{\Omega} \right)^{-1} \left( \left( \mathsf{I} - \frac{h}{4} \mathsf{\Omega} \right) - \frac{h}{2}\mathsf{\Omega} \right)
%\end{split}
%\end{equation}
%

If we approximate $\mathsf{C}^\frac12 = \frac{1}{2}(\1+\mathsf{C})$,
the half-step collisions are given by \eqref{eq:half-col-approx} and
read
%
%\begin{equation}\label{eq:half-col-approx}
%\mathsf{C}^\frac12 \fvec = \fvec^\eq + \sqrt{1-\Lambdaop} \left[ \fvec - \fvec^\eq \right] \approx \fvec^\eq + (1-\frac{1}{2}\Lambdaop) \left[ \fvec - \fvec^\eq \right] = \frac{1}{2} ( \mathsf{I} + \mathsf{C} ) \fvec ,
%\end{equation}
%%
%which leads to the approximate relation between half-step and
%full-step collisions $\Lambdaop_\frac12=\Lambdaop/2$. Application to
%$\fvec(\xvec,t)$ gives
%
\begin{equation}\label{eq:half-collision}
\begin{split}
\mathsf{C}^\frac{1}{2} \fvec(\xvec,t) 
&= \fvec(\xvec,t) - \left( \1 + \frac{h}{2} \Omegaop \right)^{-1} \frac{h}{2} \Omegaop \left[ \fvec(\xvec,t) - \fvec^{\eq}(\xvec,t) \right] + O(h^2) .\\
%&= \fvec(\xvec,t) - \frac{h}{2} \mathsf{\Lambdaop} \left[ \fvec(\xvec,t) - \fvec^{\eq}(\xvec,t) \right].
\end{split}
\end{equation}
%
In this approximation, as pointed out by Dellar \cite{Dellar2013}, the
second half-step collision \eqref{eq:trafoHe} coincides with the
transformation introduced by He et al. \cite{He1998}. This
demonstrates that the different interpretations of the operator
splitting lead to equivalent LB schemes. More precisely, the LB
variables $\tilde f_i$ in \eqref{eq:LBE} correspond in fact to the
intermediate post-streaming distributions $\bar f_i$ in the splitting
scheme, and the transformation \eqref{eq:trafoHe} is necessary to
complete the full time-step.
%
It is to be noted that the approximation in
Eq.~\eqref{eq:half-col-approx} is insufficient for global second order
accuracy if the half-step collisions are applied in every step, since
the error is $O(h^2)$ per step. Hence, the choice $\lambda_\frac12 =
\lambda/2$ does not lead to the same viscosity as the standard scheme,
which can be easily seen from Eq.~\eqref{eq:viscosity}. The difference
can be traced back to the shift $h/4$ in Eq.~\eqref{eq:BGKlambda}
which is a consequence of the splitting of the collisions into two
half-steps. The correct relaxation parameter $\lambda_\frac12$ can
instead be directly determined from Eq.~\eqref{eq:viscosity} or
\eqref{eq:BGKlambda}. It should also be noted that it is not
straightforward to invert \eqref{eq:half-collision} if the collision
operator includes fluctuations.

The force operator $\mathsf{F}$ is readily included in the splitting
scheme $\mathsf{F}^\frac12 \mathsf{C}^\frac12 \mathsf{S}
\mathsf{C}^\frac12 \mathsf{F}^\frac12$. Equations \eqref{eq:csc} for
the intermediate steps become
%
\begin{subequations}
\begin{align}
&\fvec(\xvec,t) =
\begin{multlined}[t]
%\fvec(\xvec+h\cvec,t+h) = 
\fvec^\eq(\rho,\uvec(t-\frac{h}{2})) 
+ (\1-\Lambdaop_\frac12) \bigg[ \bar \fvec(\xvec,t) 
- \fvec^\eq(\rho, \uvec(t-\frac{h}{2}) \bigg] \\
+ \frac{h}{2} \vec{G}(t-\frac{h}{4})
\end{multlined} \label{eq:force-second-halfstep} \\
&\bar f_i(\xvec+h\cvec_i,t+h) = f_i^*(\xvec,t) \label{eq:force-streaming}\\
&\fvec^*(\xvec,t) = \fvec^\eq(\rho,\uvec(t)) + (\1 - \Lambdaop_\frac12) \left[ \fvec(\xvec,t) - \fvec^\eq(\rho,\uvec(t)) \right] + \frac{h}{2} \vec{G}(t+\frac{h}{4}) ,
\label{eq:force-first-halfstep}
\end{align}
\end{subequations}
%
where we have used midpoint forces in both forcing steps. As before,
the relation to the standard LB scheme can be elucidated in terms of
the post-streaming variables $\bar f_i(\xvec+h\cvec_i,t+h)$. Inserting
\eqref{eq:force-second-halfstep} and \eqref{eq:force-first-halfstep}
into \eqref{eq:force-streaming} we obtain 
%
\begin{equation}
\begin{split}
&\bar f_i(\xvec+h\cvec,t+h) \\
&= f_i^\eq(\rho,\uvec(t)) + \sum_j (\1 - \Lambdaop_\frac12)_{ij} \left[ f_j(\xvec,t) - f_j^\eq(\rho,\uvec(t)) \right] + \frac{h}{2} G_i(t+\frac{h}{4})\\
%&= \fvec^\eq(\rho,\uvec(t)) + (1 - \Lambdaop_\frac12) \left[ \fvec(\xvec,t) - \fvec^\eq(\rho,\uvec(t-\frac{h}{2})) - \frac{h}{2} \vec{G}(t+\frac{h}{4}) \right] + \frac{h}{2} \vec{G}(t+\frac{h}{4})\\
&= f_i^\eq(\rho,\uvec(t)) + \sum_j (\1 - \Lambdaop_\frac12)^2_{ij} \left[ \bar f_j(\xvec,t) - f_j^\eq(\rho,\uvec(t-\frac{h}{2})) \right] + \frac{h}{2} G_i(t+\frac{h}{4})\\
&= 
\begin{multlined}[t]
f_i^\eq(\rho,\uvec(t)) + \sum_j (\1 - \Lambdaop_\frac12)^2_{ij} \left[ \bar f_j(\xvec,t) - f_j^\eq(\rho,\uvec(t)) \right] \\
+ \sum_j (\1-\Lambdaop_\frac12)^2_{ij} \frac{h}{2} G_j(t-\frac{h}{4}) 
+ \frac{h}{2} G_i(t+\frac{h}{4}) .
\end{multlined}
\end{split}
\end{equation}
%
Assuming that the force term $G_i$ can be expanded around $t$, we find
up to terms of $O(h^3)$
%
\begin{equation}
\begin{split}
&\bar f_i(\xvec+h\cvec_i,t+h) \\
&= f_i^\eq(\rho,\uvec(t)) + 
\begin{multlined}[t]
\sum_j (\1 - \Lambdaop_\frac12)^2_{ij} \left[ \bar f_j(\xvec,t) - f_j^\eq(\rho,\uvec(t)) \right] \\
+ \sum_j \left( \1 + (\1-\Lambdaop_\frac12)^2 \right)_{ij} \frac{h}{2} G_j(t) + O(h^3)
\end{multlined} \\
&= f_i^\eq(\rho,\uvec(t)) + 
\begin{multlined}[t]
\sum_j (\1 - \Lambdaop)_{ij} \left[ \bar f_j(\xvec,t) - f_j^\eq(\rho,\uvec(t)) \right] \\
+ \sum_j (\1 - \frac{1}{2}\Lambdaop)_{ij} h G_j(t) + O(h^3) .
\end{multlined} \\
\end{split}
\label{eq:force-transformation}
\end{equation}
%
Equation \eqref{eq:force-transformation} contains the familiar
transformation of the force term $\vec{G}$ that is used in standard LB
schemes to maintain second-order accuracy
\cite{He1998a,Ladd2001,Guo2002,Duenweg2008}.  In the splitting scheme
it arises naturally without the need to know the transformation
a-priori.

Another aspect that is illuminated by the derivation is the
redefinition of the hydrodynamic velocity in the presence of external
forces \cite{Ladd2001,Duenweg2008}
%
\begin{equation}
\rho \uvec = \sum_i \bar f_i \cvec_i \cvec_i + \frac{h}{2} \Fvec ,
\end{equation}
%
which is simply a consequence of the second application of the force
operator $\mathsf{F}^\frac12$ after the streaming step.

Although here the primary purpose of the collide-stream-collide scheme
is to elucidate the operator splitting approach, we remark that it is
straightforward to implement using a combination of the familiar
``push'' and ``pull'' schemes \cite{Wittmann2013}. A simple
pseudo-code is listed in Algorithm~\ref{alg:csc}. If one uses a
memory-efficient data layout for intermediate results, there is only a
small overhead in form of the floating point operations for the second
relaxation step. This may be almost negligible because the performance
of the LBM on modern computing architectures is mainly limited by
memory bandwidth, and the bottleneck regarding lattice site updates
per second is really the streaming step.  Nevertheless, it may seem
that in terms of floating point operations, the collide-stream-collide
LB scheme is less efficient due to the extra collision step. On the
other hand, the standard LB scheme requires extra effort to evaluate
the viscous stress, cf. Eq.~\eqref{eq:viscstress}. As long as the
stress is only needed at the end of the calculation, the standard LB
scheme may be preferable. However, if a coupling scheme involves the
fluid stress, its computation is needed after every update. In this
case the collide-stream-collide is almost equally efficient, and it
has the additional advantage that the viscous stress is directly
available as the local non-equilibrium stress.
%
%\comment{A detailed investigation of the accuracy and efficiency of
%collide-stream-collide lattice Boltzmann is left for future work.}



\section{Coupling methods for hybrid fluid-particle systems}

We now turn to methods for coupling immersed particles to the lattice
Boltzmann equation that allow to simulate complex fluids. Here, we
focus on objects with internal degrees of freedom, such as polymers or
vesicles. Unlike colloids, they can not easily be coupled to the fluid
by simple boundary conditions due to their complex and changing
shape. The most common methods employed in such cases are the immersed
boundary method (IBM), force coupling (FC), and external boundary
force (EBF).



\subsection{Immersed boundary method (IBM)}
\label{sec:IBM}

In the immersed boundary method \cite{Peskin2002}, a fluid-object
interface is represented by a set of Lagrangian nodes with positions
$\Xvec_i$ that are advected with the flow. Interactions on the
Lagrangian grid are applied as body forces to the fluid. The positions
of the Lagrangian nodes are additional degrees of freedom that are
updated according to
%
\begin{equation}
\frac{d}{d t} \Xvec_i(t) = \mathcal{P} \Xvec_i(t) = \uvec(\Xvec_i,t) ,
\end{equation}
%
where $\mathcal{P}$ is the propagation operator and $\uvec(\Xvec_i,t)$
is the flow velocity at the Lagrangian position $\Xvec_i$. 
%
The latter requires a mapping between Lagrangian positions and
Eulerian grid coordinates
%
\begin{equation}\label{eq:IBM-interpolation}
\uvec(\Xvec_i,t) = \Iop_a[\Xvec_i(t)] \, \uvec (\Xvec_i) = \sum_\xvec \uvec(\xvec,t) \delta_a(\Xvec_i(t)-\xvec) ,
\end{equation}
%
where $\delta_a$ is a discrete Dirac delta function on a
lattice with lattice spacing $a$, and $\Iop_a[\Xvec_i]$ is an
interpolation operator based on the Lagrangian positions $\Xvec_i$.
%
Conversely, the forces arising between the Lagrangian nodes are spread
to the grid points with the adjoint operator $\Iop_a^*$
%
\begin{equation}\label{eq:IBM-spreading}
\Fvec(\xvec,t) = \Iop_a^*[\Xvec_i(t)] \, \Fvec_i (\xvec) = \frac{1}{a^3} \sum_{i} \Fvec_i(t) \delta_a(\Xvec_i(t)-\xvec) .
\end{equation}
%
Note that $\Fvec(\xvec,t)$ is a force density while $\Fvec_i(t)$ is a
force, hence the factor $a^{-3}$ on the right hand side of
\eqref{eq:IBM-spreading}.
%
The construction of the discrete delta function $\delta_a$ is
described in detail in \cite{Peskin2002} and is briefly summarized in
\ref{sec:delta}.

The evolution of the complete system is described by
%
\begin{equation}
\frac{\d}{\d t} \{ \fvec(\xvec,t), \Xvec_i(t) \} = \left[ \mathcal{S} + \mathcal{C} + \mathcal{F} + \mathcal{P} \right] \{ \fvec(\xvec,t), \Xvec_i(t) \}.
\end{equation}
%
Since the interpolation operator $\Iop_a[\Xvec_i(t)]$ depends
on the positions of the Lagrangian nodes, $\mathcal{P}$ does not
commute with any of the operators $\mathcal{S}$, $\mathcal{C}$ and
$\mathcal{F}$.
%
Therefore, the commonly used IBM algorithm
%
\begin{equation}
\{ \fvec(\xvec,t+h), \Xvec_i(t+h) \} = \left[ \mathsf{P}\mathsf{S}\mathsf{C}\mathsf{F} \right] \{ \fvec(\xvec,t), \Xvec_i(t) \}
\end{equation}
%$\mathsf{P}\mathsf{S}\mathsf{C}\mathsf{F}$ 
is only first-order accurate. It is straightforward to write down a
formally second order scheme:
%
\begin{equation}\label{eq:IBM-scheme}
\{ \fvec(\xvec,t+h), \Xvec_i(t+h) \} = \left[ \mathsf{P}^\frac12 \mathsf{F}^\frac12 (\mathsf{C}^\frac12 \mathsf{S} \mathsf{C}^\frac12) \mathsf{F}^\frac12 \mathsf{P}^\frac12 \right] \{ \fvec(\xvec,t), \Xvec_i(t) \} .
\end{equation}
%
During the position updates $\mathsf{P}^\frac12$, the velocity is
constant but the interpolation $\Iop_a[\Xvec_i(t)]$ changes. Therefore
a simple forward Euler update is insufficient, and the Crank-Nicolson
rule becomes implicit. While one possibility is to use the midpoint
rule in each of the two half-steps $\mathsf{P}^\frac12$, we suggest
here a simpler explicit scheme that uses a predictor step for the
positions after the first half-step. Then we use a backward Euler in
the first position half-step and a forward Euler in the second half
step:
%
\begin{equation}
\begin{split}
\tilde \Xvec_i(t+\frac{h}{2}) &= \Xvec_i(t) + \frac{h}{2} \Iop_a[\Xvec_i(t)]\uvec(t) , \\
\Xvec_i(t+\frac{h}{2}) &= \Xvec_i(t) + \frac{h}{2} \Iop_a[\tilde\Xvec_i(t+\frac{h}{2})]\uvec(t) , \\
\Xvec_i(t+h) &= \Xvec_i(t+\frac{h}{2}) + \frac{h}{2} \Iop_a[\tilde\Xvec_i(t+\frac{h}{2})]\uvec(t+h) .
\end{split}
\end{equation}
%
We get an overall position update
%
\begin{equation}
\Xvec_i(t+h) = \Xvec_i(t) + h \Iop_a[\tilde\Xvec_i(t+\frac{h}{2})] \frac{\uvec(t) + \uvec(t+h)}{2} ,
\end{equation}
%
which is a trapezoidal quadrature of Eq.~\eqref{eq:IBM-interpolation}.
%and is second order time-accurate despite the lower order of the
%sub-steps. \comment{Is this true?}  
This scheme can be seen as the adaption of the formally second-order
method in Refs.~\cite{Lai2000,Mori2008} to the lattice Boltzmann
method. If the body force $\Fvec(\xvec,t)$ depends on the positions
$\Xvec_i$ but not on the flow velocity $\uvec$, the body forces in the
first and second half-steps $\mathsf{F}^\frac12$ are the same. The
application of the body forces thus relates to a midpoint scheme which
is interleaved with the LB collisions and streaming.
%
\begin{equation}
\begin{split}
\rho \uvec(t+h) &= \rho \uvec(t+\frac{h}{2}) 
+ \frac{h}{2} \Iop_a^*[\tilde\Xvec_i(t+\frac{h}{2})] \, \Fvec_i(t+\frac{h}{2}) \\
&
\begin{multlined}[t]
= \left[ \mathsf{C}^\frac12 \mathsf{S} \mathsf{C}^\frac12 \right]
\left( \rho \uvec(t) + \frac{h}{2}  \Iop_a^*[\tilde\Xvec_i(t+\frac{h}{2})] \, \Fvec_i
(t+\frac{h}{2}) \right) \\
+ \frac{h}{2} \Iop_a^*[\tilde\Xvec_i(t+\frac{h}{2})] \, \Fvec_i(t+\frac{h}{2}) . 
\end{multlined}
\end{split}
\end{equation}
%
The scheme in Eq.~\eqref{eq:IBM-scheme} is easily implemented by
combining a Verlet-like update for the Lagrangian nodes $\Xvec_i$ with
any second order accurate LB scheme. A pseudo-code is listed in
Algorithm \ref{alg:ibm}.



\subsection{Force coupling}
\label{sec:force-coupling}

The immersed boundary method works well if the dynamics of the fluid
and the solute objects evolve on similar time scales. In soft matter
systems, the internal motion of the solutes can be much faster than
the fluid motion and hence needs to be resolved separately. This can
be done by coupling LB to a molecular dynamics simulation where the
immersed particles evolve according to Newton's equations of motion
\cite{Ahlrichs1999}. The particle positions $\rvec_i$ are now
Lagrangian points that move with their own velocity $\vvec_i$, and the
intrinsic interactions $\Fvec_i$ enter Newton's equation of motion for
the particles.

The fluid-particle interactions are incorporated by a coupling force that
constrains the fluid flow to the motion of particles and vice-versa.
%
This coupling force is modeled as a drag force proportional to the
relative velocity of the particles with respect to the flow
\cite{Ahlrichs1999}
%
\begin{equation}\label{eq:viscous-coupling}
\begin{split}
\Fvec^{h}_i(t) &= - \gamma \left[ \vvec_i(t) - \uvec(\rvec_i,t) \right] + \zetavec_i ,
%= - \gamma \left[ \vvec_i(t) - (\Iop_a[\rvec_i(t)] \, \uvec)(\rvec_i,t) \right] + \zetavec_i , \\
%&= - \gamma \left [ \vvec_i(t) - \sum_\xvec \uvec(\xvec,t) \, a^3 \delta_a(\rvec_i(t)-\xvec) \right] + \zetavec_i 
\end{split}
\end{equation}
%
where $\gamma$ is a friction constant, $\zetavec_i$ is a random force
required to satisfy the fluctuation-dissipation theorem, and
$\uvec(\rvec_i,t)=\Iop_a[\rvec_i(t)]\uvec(\rvec_i,t)$ is the fluid
velocity at the particle's position.
%
Since the coupling force must not change the total momentum of the
combined fluid-particle system, an opposite force needs to be applied
to the fluid in terms of body forces on the lattice sites within the
interpolation range around the particle
%
\begin{equation}
\begin{split}
\Fvec^{h}(\xvec,t) &= - \Iop_a^*[\rvec_i(t)] \Fvec_i^h(t) . %\Iop_a^*[\rvec_i(t)] \, \left( \gamma \left[ \vvec_i(t) - (\Iop_a[\rvec_i(t)] \, \uvec)(\rvec_i,t) \right] - \zetavec_i \right) (\xvec,t)  . \\
%&= -\gamma (\Iop_a^*[\rvec_i(t)] \Iop_a[\rvec_i(t)] \uvec) (\xvec,t) - \gamma \uvec(\xvec,t) + \Iop_a^*[\rvec_i(t)] \, \left( \gamma \vvec_i - \zetavec_i \right) (\xvec,t)\\
%&= \sum_i \delta_a(\rvec_i(t)-\xvec) \left[ \gamma \vvec_i(t) - \gamma \sum_{\xvec'} \uvec(\xvec',t) \, a^3 \delta_a(\rvec_i-\xvec') + \zetavec_i \right]
\end{split}
\end{equation}
%
The interpolation and spreading operators are the same as in the
immersed boundary method. The equations of motion for the particles are
%
\begin{subequations}
\begin{align}
\frac{d}{d t}\rvec_i(t) &= \mathcal{P} \rvec_i(t) = \vvec_i(t) , \\
\frac{d}{d t}\vvec_i(t) &= \mathcal{F} \vvec_i(t) = \frac{1}{m_i} \left( \Fvec^{h}_i + \Fvec_i^\text{int}  \right) ,
\end{align}
\end{subequations}
%
where $\Fvec_i^\text{int}$ are intrinsic interactions between the
particles.

The evolution of the coupled system is now given by
%
\begin{equation}
\frac{\d}{\d t} \left\{ \fvec(\xvec,t), \rvec_i(t), \vvec_i(t) \right\} = \left[ \mathcal{S} + \mathcal{C} + \mathcal{F}  + \mathcal{P} \right] \left\{ \fvec(\xvec,t), \rvec_i(t), \vvec_i(t) \right\} ,
\end{equation}
%
where the overall force operator $\mathcal{F}$ acts simultaneously on
the particle velocities and on the fluid flow,
cf.~Eq.~\eqref{eq:LB-force-moments},
%
\begin{subequations}
\begin{align}
&\mathcal{F} \vvec_i(t)     = \Fvec_i^h(t) , \\ 
&\mathcal{F} \uvec(\xvec,t) = - \Iop_a^*[\rvec_i(t)] \Fvec_i^h(t) .
\end{align}
\end{subequations}
%
As before, we can apply the Trotter-Suzuki decomposition formula to
obtain a splitting scheme \cite{Suzuki1985}
%
\begin{equation}\label{eq:fc-splitting}
\left\{ \fvec, \rvec_i, \vvec_i \right\}(t+h) = \Big[ \mathsf{P}^\frac12 \mathsf{F}^\frac12 ( \mathsf{C}^\frac12 \mathsf{S} \mathsf{C}^\frac12 ) \mathsf{F}^\frac12 \mathsf{P}^\frac12 \Big] \left\{ \fvec, \rvec_i, \vvec_i \right\}(t) .
\end{equation}
This is basically a Verlet algorithm with a lattice Boltzmann update
in the middle of the velocity update.
%
Since the hydrodynamic coupling forces are dependent on both the
particle velocities and the fluid flow, the equations of motion are
coupled
%
\begin{subequations}\label{eq:fc-vel}
\begin{align}
\frac{\d}{\d t} \vvec_i(t) &= - \frac{1}{m_i} \left[ \gamma \left( \vvec_i(t) - \uvec(\rvec_i,t) \right) - \zetavec_i - \Fvec_i^\text{int} \right] , \\%\label{eq:fc-vel}\\
\frac{\d}{\d t} \uvec(\rvec_i,t) &= \frac{1}{\rho a^3} \left[ \gamma \left( \vvec_i(t) - \uvec(\rvec_i,t) \right) - \zetavec_i \right] .
\end{align}
\end{subequations}
%
An accurate midpoint scheme for the velocities and the fluid flow is
%
\begin{subequations}\label{eq:fc-midpoint} 
\begin{align}
\vvec_i(t+h) &= \vvec_i(t) +\frac{h}{m_i} \Fvec_i^h(t+\frac{h}{2}) + \frac{h}{m_i} \Fvec_i^\text{int} , \\
%- \frac{\gamma h}{m_i} \left[ \tilde\vvec_i(t+\frac{h}{2}) - \Iop_a[\rvec_i] \tilde\uvec(t+\frac{h}{2}) -\frac1\gamma \zetavec_i \right] + \frac{h}{m_i} \Fvec_i^\text{int} , \\
\uvec(\rvec_i,t+h) &= \uvec(\rvec_i,t) - \frac{h}{\rho} \Fvec^h_i(t+\frac{h}{2})  . %\label{eq:fc-midpoint2} ,
%+ \frac{\gamma h}{m_i} \Iop_a^*[\rvec_i] \left[ \tilde\vvec_i(t+\frac{h}{2}) - \Iop_a[\rvec_i] \tilde\uvec(t+\frac{h}{2}) -\frac1\gamma \zetavec_i \right] 
\end{align}
\end{subequations}
%
%where the midpoint force $\Fvec_i^h(t+\frac{h}{2})$ could be
%calculated using predictor steps for the velocities.
%
In Ref. \cite{Ladd2009}, a constant flow velocity $\uvec(\rvec_i,t)$
is assumed to determine the viscous coupling force. This is justified
for slowly varying flow fields that are only slightly perturbed by the
coupling, but strictly speaking it does not give an $O(h^2)$
approximation of Eq.~\eqref{eq:fc-vel}.
%
%The result in \cite{Ladd2009} approximates the exact solution of
%Eq.~\eqref{eq:fc-vel} up to and including terms of $O(h^2)$
%%
%%\begin{equation}
%%\d_t \vvec_i = - \frac{\gamma}{m_i} \left[ \vvec_i(t) - \Iop_a[\rvec_i] \uvec - \frac{1}{\gamma} ( \zetavec_i + \Fvec_i^\text{int} ) \right]
%%\end{equation}
%\begin{equation}
%\vvec_i(t+h) = \vvec_i(t) e^{-\alpha} + \left( 1 -e^{-\alpha} \right)
%\left( \uvec(\rvec_i,t) + \frac{1}{\gamma} ( \zetavec_i +
%\Fvec_i^\text{int} ) \right).
%\end{equation}
%
%Eq.~\eqref{eq:fc-midpoint3} is up to $O(h^2)$ equivalent to the
%approximation suggested in Ref.~\cite{Ladd2009}, where the intrinsic
%forces are applied separately. This corresponds to another splitting
%of the force operator and is compatible with second-order
%time-accuracy.
%
Therefore, we directly apply Crank-Nicolson discretization here
%
%Here, we assume that the flow velocity $\tilde\uvec(t+\frac{h}{2})$ is
%constant during the hydrodynamic interaction \cite{Ladd2009}.
%
%However, these are coupled equations and it is not immediately obvious
%how to determine $\vvec_i(t+\frac{h}{2})$ and $\uvec(t+\frac{h}{2})$
%such that exact overall momentum conservation and the FDT can be
%satisfied.
%
%If we assume that the momentum exchange due to the coupling force
%takes place instantaneously at the half-step, we get the implicit
%equation
%%
%\begin{align}
%%\vvec_i(t+\frac{h}{4}) &= \vvec_i(t) + \frac{h}{4m_i} \Fvec_i^\text{int}\\
%\tilde \vvec_i(t+\frac{h}{2}) %&= \vvec_i(t+\frac{h}{4}) + \frac{h}{2m_i} \Fvec_i^h(\tilde\vvec_i(t+\frac{h}{4})) \\
%&=  \vvec_i(t) + \frac{h}{2m_i} \Fvec_i^\text{int} - \frac{\gamma h}{m_i} \left[ \tilde\vvec_i(t+\frac{h}{2}) - \Iop_a[\rvec_i] \tilde\uvec(t+\frac{h}{2}) - \frac1\gamma \zetavec_i \right] \\
%%\vvec_i(t+\frac{h}{2}) &= \tilde\vvec_i(t+\frac{h}{4}) + \frac{h}{4m_i} \Fvec_i^\text{int}
%\end{align}
%
%We can then apply Crank-Nicolson discretization
%
\begin{subequations}\label{eq:fc-cn}
\begin{align}
\begin{multlined}[t]
\frac{\vvec_i(t+h) - \vvec_i(t)}{h} = 
- \frac{\gamma}{m_i} \frac{\vvec_i(t+h) - \uvec(\rvec_i,t+h)}{2} \\
- \frac{\gamma}{m_i} \frac{\vvec_i(t) - \uvec(\rvec_i,t)}{2}
+ \frac{1}{m_i} \left( \Fvec_i^\text{int} +\zetavec_i \right) ,
\end{multlined} \\
\begin{multlined}[t]
\frac{\uvec(\rvec_i,t+h) - \uvec(\rvec_i,t)}{h} = 
\frac{\gamma}{\rho a^3} \frac{\vvec_i(t+h) - \uvec(\rvec_i,t+h)}{2} \\
+ \frac{\gamma}{\rho a^3} \frac{\vvec_i(t) - \uvec(\rvec_i,t)}{2}
- \frac{1}{\rho a^3} \zetavec_i ,
\end{multlined}
\end{align}
\end{subequations}
%
%\begin{align}\label{eq:fc-cn}
%\vvec_i(t+h) = \vvec_i(t) - \frac{\gamma h}{m_i} \frac{\vvec_i(t+h)+\vvec_i(t)}{2} + \frac{h}{m_i} \left( \gamma \Iop_a[\rvec_i] \uvec + \Fvec_i^\text{int} + \zetavec_i \right) ,
%\end{align}
%
%\begin{equation}
%\vvec_i(t+h) \left( 1 + \frac{\gamma h}{2m_i} \right) = \vvec_i(t) \left( 1 - \frac{\gamma h}{2m_i} \right) + \frac{h}{m_i} \left( \Fvec_i^\text{int} + \zetavec_i + \gamma \Iop_a[\rvec_i] \uvec \right) 
%\end{equation}
%
and obtain
%
\begin{align}
\vvec_i(t+h) &= \vvec_i(t) - \frac{\gamma^{-1}\alpha}{1+\frac{\alpha}{2}+\frac{\beta}{2}} \left[ \gamma \left( \vvec_i(t) - \uvec(\rvec_i,t) \right) - \zetavec_i - \frac{2+\beta}{2} \Fvec_i^\text{int} \right] , \\
\uvec(\rvec_i,t+h) &= \uvec_i(\rvec_i,t) + \frac{\gamma^{-1}\beta}{1+\frac{\alpha}{2}+\frac{\beta}{2}} \left[ \gamma \left( \vvec_i(t) - \uvec(\rvec_i,t) \right) - \zetavec_i + \frac{\alpha}{2} \Fvec_i^\text{int} \right] .
\end{align}
%
%\begin{equation}\label{eq:fc-midpoint3}
%\begin{split}
%\vvec_i(t+h) &= \vvec_i(t) - \frac{\alpha}{1 + \frac{\alpha}{2}} \left( \vvec_i(t) - \Iop_a[\rvec_i] \uvec - \frac{1}{\gamma} ( \Fvec_i^\text{int} + \zetavec_i ) \right) \\
%&= \vvec_i(t) + \frac{\alpha}{1 + \frac{\alpha}{2}} \frac{ \Fvec_i^h(t) + \Fvec_i^\text{int} }{\gamma} ,
%\end{split}
%\end{equation}
%
where the parameters $\alpha=\frac{\gamma h}{m_i}$ and
$\beta=\frac{\gamma h}{\rho a^3}$ were introduced.
%, and $\Fvec_i^h(t)$ is the coupling force \eqref{eq:viscous-coupling}
%evaluated with the velocities $\vvec_i(t)$ and $\uvec(\rvec_i,t)$.
%
Requiring consistency between \eqref{eq:fc-midpoint} and
\eqref{eq:fc-cn} we obtain the midpoint force
%
%\begin{equation}\label{eq:modified-force}
%\begin{split}
%\Fvec_i^h(t+\frac{h}{2}) &= \frac{m_i}{h} \left[ \vvec_i(t+h) - \vvec_i(t) \right] - \Fvec_i^\text{int} \\
%%&= \frac{1}{1+\frac{\alpha}{2}} \left( \Fvec_i^h(t) + \Fvec_i^\text{int} \right) - \Fvec_i^\text{int} \\
%%&= \frac{1}{1+\frac{\alpha}{2}} \Fvec_i^h(t) - \frac{\frac{\alpha}{2}}{1+\frac{\alpha}{2}} \Fvec_i^\text{int} \\
%&= \frac{1}{1+\frac{\alpha}{2}} \left( \Fvec_i^h(t) -\frac{\alpha}{2} \Fvec_i^\text{int} \right) .
%\end{split}
%\end{equation}
%
\begin{equation}\label{eq:modified-force}
\begin{split}
\Fvec_i^h(t+\frac{h}{2}) &= - \frac{1}{1+\frac{\alpha}{2}+\frac{\beta}{2}} \left[ \gamma \left( \vvec_i(t) - \uvec(\rvec_i,t) \right) - \zetavec_i + \frac{\alpha}{2} \Fvec_i^\text{int} \right] .
\end{split}
\end{equation}
%
The denominator $1+\frac{\alpha}{2}+\frac{\beta}{2}$ has the typical
form of a discrete correction, cf. Eq.~\eqref{eq:collision-operator},
and takes into account that both $\vvec_i(t)$ and $\uvec(\rvec_i,t)$
change during the time step.
%
%\comment{Can this be interpreted somehow else?}
%
%\begin{equation}
%\vvec_i(t+\frac{h}{2}) = \frac{m_i}{\gamma h} \left[ \vvec_i(t) ( 1 - e^{-\frac{\gamma h}{m_i}} ) + \left( \frac{\gamma h}{m_i} - 1 + e^{-\frac{\gamma h}{m_i}}\right) \left( \Iop_a[\rvec_i] \tilde\uvec(t+\frac{h}{2}) + \frac{1}{\gamma} ( \zetavec_i +\Fvec_i^\text{int} ) \right) \right]
%\end{equation}
%
%Eq.~\eqref{eq:fc-midpoint3} leads to more accurate integration of the
%coupling forces \cite{Ladd2009} and can also be interpreted as a
%forward Euler step with the force
%\begin{equation}
%\Fvec_i = \Fvec_i^\text{int} + \Fvec_i^h(t+\frac{h}{2}) = \frac{1}{1+\frac{\alpha}{2}} \left(
%\Fvec_i^\text{int} + \Fvec_i^h(t) \right).
%\end{equation}
Eq.~\eqref{eq:modified-force} is straightforward to implement and
can lead to more accurate integration of the coupling forces
\cite{Ladd2009}.
%
%To ensure that the coupling force does not change the total momentum
%of the system, its negative is spread to the fluid 
%%
%\begin{align}\label{eq:fluid-force}
%%\Fvec_i^h(t+\frac{h}{2}) &= - \gamma \left[ \tilde\vvec(t+\frac{h}{2}) - \Iop_a[\rvec_i]\tilde\uvec(t+\frac{h}{2}) \right] + \zetavec_i \\
%\Fvec^h(\xvec,t+\frac{h}{2}) &= - \Iop_a^*[\rvec_i] \Fvec_i^h(t+\frac{h}{2})(\xvec) .
%\end{align}
%%
%This body force is the action of the operator $\mathcal{F}$ on the LB
%variables $f_i(\xvec,t+h)=f_i(\xvec,t) + G_i(t+h/2)$ with the moments
%of $G_i$ given by Eq.~\eqref{eq:force-moments}.
%
%Eq.~\eqref{eq:fc-midpoint}
%and \eqref{eq:fc-midpoint2} together with the midpoint coupling force
%\eqref{eq:modified-force} specify the action of the force operator
%$\mathsf{F}$.  
%
Note that the formulas
\eqref{eq:fc-midpoint}-\eqref{eq:modified-force} refer to the full
time step operator $\mathsf{F}$, and a half-step $\mathsf{F}^\frac12$
is achieved by substituting $h$ by $\frac{h}{2}$. For a full time step
in the splitting scheme \eqref{eq:fc-splitting}, the procedure is
applied before and after the LB update. A pseudo-code for force
coupling is listed in Algorithm~\ref{alg:fc}.

It is worthwhile to note that the viscous force coupling does not lead
to an instantaneous stick boundary condition. To see this, we
calculate the difference between the particle velocity and fluid flow
after the force application
%
\begin{multline}
\vvec_i(t+h) - \uvec(\rvec_i,t+h) = \frac{1 - \frac{\alpha}{2} -
\frac{\beta}{2}}{1+\frac{\alpha}{2}+\frac{\beta}{2}} \left( \vvec_i(t)
- \uvec(\rvec_i,t) \right) \\
+ \frac{\alpha+\beta}{1+\frac{\alpha}{2}+\frac{\beta}{2}}
\frac{\zetavec_i}{\gamma} +
\frac{\alpha}{1+\frac{\alpha}{2}+\frac{\beta}{2}} 
\frac{\Fvec_i^\text{int}}{\gamma} .
\end{multline}
%
%\begin{equation}\label{eq:visc-slip}
%\begin{split}
%&\vvec_i(t+h) - \uvec(\rvec_i,t+h) \\
%&= 
%%\begin{multlined}[t]
%\vvec_i(t) - \uvec(\rvec_i,t) + \frac{\frac{h}{m_i}}{1+\frac{\alpha}{2}} \left( \Fvec_i^h(t) + \Fvec_i^\text{int} \right) + \frac{\frac{h}{\rho a^3}}{1+\frac{\alpha}{2}} \left( \Fvec_i^h(t) - \frac{\alpha}{2} \Fvec_i^\text{int} \right) \\
%%\end{multlined} \\ 
%%&= \vvec_i(t) - \uvec(\rvec_i,t) + \frac{\frac{h}{m_i} ( 1 + \frac{m_i}{\rho a^3} )}{1+\frac{\alpha}{2}} \Fvec_i^h(t) + \frac{\frac{h}{m_i} ( 1 - \frac{m_i}{\rho a^3}\frac{\alpha}{2} ) }{1+\frac{\alpha}{2}} \Fvec_i^\text{int}\\
%&= 
%%\begin{multlined}[t]
%\shoveleft{\frac{1 - \frac{\alpha}{2} ( 1 + \frac{2m_i}{\rho a^3} )}{1+\frac{\alpha}{2}} \left( \vvec_i(t) - \uvec(\rvec_i,t) \right) + \frac{\frac{h}{m_i} ( 1 - \frac{m_i}{\rho a^3}\frac{\alpha}{2} ) }{1+\frac{\alpha}{2}} \Fvec_i^\text{int}.}
%%\end{multlined} \\
%\end{split}
%\end{equation}
%
The second term stems from the thermal fluctuations and is zero on
average. 
%The condition to make both coefficients in Eq.~\eqref{eq:visc-slip}
%equal to zero is $1+\frac{m_i}{\rho a^3}=0$ which does not have a
%physical solution.
Requiring that the coefficient of the velocity difference vanishes
yields $\alpha+\beta = 2$ and an expression for the friction
%
\begin{align}\label{eq:gamma-cond}
\begin{split}
%\alpha + \beta = 2 && \Rightarrow && 
\gamma = \frac{\rho a^3}{h} \frac{2 \beta}{\alpha + \beta} = \frac{\rho a^3}{h} \frac{2}{1 + \frac{\rho a^3}{m_i}} .
\end{split}
\end{align}
%%
%\begin{align}
%\gamma = \frac{\rho a^3}{h} \frac{1}{1+\frac{\rho a^3}{2m_i}} .%&& \alpha = \frac{1}{\frac{m_i}{\rho a^3} + \frac{1}{2}} .
%\end{align}
%%
In this case, the velocity difference after force application is
%
\begin{equation}\label{eq:slip}
\vvec_i(t+h) - \uvec(\rvec_i,t+h) = \frac{h}{2m_i} \Fvec_i^\text{int} + \frac{\zetavec_i}{\gamma} ,
\end{equation}
%
which does not vanish on average.
%
The reason for this spurious velocity slip lies in the nature of the
discretization: the velocities are effectively evaluated at the
midpoint and therefore the intrinsic forces acting during the second
half-interval are not accounted for. 
%
%In the absence of intrinsic forces, the no-slip condition can be
%satisfied by the appropriate choice of the friction~$\gamma$. In the
%presence of intrinsic forces, however, the force coupling method can
%only approximately satisfy the no-slip boundary condition.
%
This does not cause a problem if one considers point-like particles
that are dragged through the fluid. In this case, one can tune
$\gamma$ to obtain a defined physical friction and an effective
particle size. The effective physical friction has to be re-normalized
due to the spatial interpolation on the discrete lattice and leads to
an effective size $R_\text{eff}$ of the
particles~\cite{Ahlrichs1999,Duenweg2008}
%
\begin{equation}
\begin{split}
\frac{1}{R_\text{eff}} &= \frac{6\pi\eta}{\gamma} + \frac{1}{g a} .
\end{split}
\end{equation}
%
%\begin{equation}
%\begin{split}
%\frac{1}{\gamma_\text{eff}} &= \frac{1}{\gamma} + \frac{1}{6\pi\eta g a} .
%%= \frac{h}{\rho a^3} \frac{1 + \frac{\rho a^3}{m_i}}{2} + \frac{1}{6\pi\eta g a}
%\end{split}
%\end{equation}
%
%We remark that the friction $\gamma$ is only a bare input
%parameter. The effective friction has to be re-normalized due to the
%spatial interpolation on the discrete lattice
%\cite{Ahlrichs1999,Duenweg2008}
%%
%\begin{equation}
%\begin{split}
%\frac{1}{\gamma_\text{eff}} &= \frac{1}{\gamma} + \frac{1}{6\pi\eta g a} .
%%= \frac{h}{\rho a^3} \frac{1 + \frac{\rho a^3}{m_i}}{2} + \frac{1}{6\pi\eta g a}
%\end{split}
%\end{equation}
%
%\begin{equation}
%\begin{split}
%\frac{1}{R_\text{eff}} &= \frac{6 \pi \eta}{\gamma} + \frac{1}{ga}
%= \frac{6 \pi \eta h}{\rho a^3} \frac{1 + \frac{\rho a^3}{m_i}}{2} + \frac{1}{ga}
%\end{split}
%\end{equation}
%
Moreover, since the model is based on point particles, extended
objects will necessarily be under-resolved when represented by a
single particle. Therefore it is necessary to add fluctuations
$\zetavec_i$ to the coupling force in order to satisfy the fluctuation
dissipation relation \cite{Ahlrichs1999,Duenweg2008}.

In the case of more complicated objects like cells and vesicles with
specific shapes, it may me more appropriate to represent the object by
a set of surface points that are coupled to the fluid. The spurious
slip velocity \eqref{eq:slip} may then be undesirable. A remedy for
this situation is a modification of the external boundary force as
described in the following section.



\subsection{External boundary force (EBF)}
\label{sec:ebf}

The modeling of the fluid-particle interactions as a viscous force has
proved successful for simulating polymer chains in solution, where the
individual beads of the chain experience the viscous
coupling~\cite{Ahlrichs1999,Pham2009,Ladd2009}, but it may be less
appropriate for extended objects. In this case, the fluid velocity at
the surface is typically required to satisfy a no-slip boundary
condition, i.e., the fluid and the object should move with the same
velocity. Therefore, the fluid-particle interaction can also be
modeled as an external boundary force that forces the fluid to move
with the surface velocity \cite{Goldstein1993,Wu2010}. In the latter
references, the interaction force acting at a surface point is
%
\begin{equation}\label{eq:ebf}
\Fvec^h_i(\rvec_i,t) = - \frac{\rho a^3}{h} \left( \vvec_i(t) - \uvec(\rvec_i,t) \right) ,
\end{equation}
%
which is equivalent to the force coupling \eqref{eq:viscous-coupling}
with $\gamma=\frac{\rho a^3}{h}$. In this section we skip the
stochastic force $\zetavec_i$ which is straightforward to include if
demanded by the fluctuation dissipation relation for the immersed
object.

As we have seen above, however, the force in \eqref{eq:ebf} can not
lead to the no-slip boundary condition if intrinsic forces are
present.
%The reason is that the fluid is constrained to the
%surface velocity, but not vice-versa. While this may be appropriate
%for stationary rigid walls, it is less suited for moving or deformable
%objects such as vesicles and cells. 
The reason is that, in contrast to the immersed boundary method, the
intrinsic forces are applied to the particles and accelerate them
relative to the fluid. In order to enforce the no-slip boundary
condition even in the presence of intrinsic forces, we require
%
\begin{equation}
\begin{split}
\vvec_i(t+h) &= \uvec(\rvec_i,t+h) , \\
\vvec_i(t) + \frac{h}{m_i} \left[ \Fvec_i^h(t+\frac{h}{2}) + \Fvec_i^\text{int} \right]
&= \uvec(\rvec_i,t) - \frac{h}{\rho a^3} \Fvec_i^h(t+\frac{h}{2}) .
%= \uvec(\rvec_i,t+h)
\end{split}
\end{equation}
%
Solving for the fluid-particle interaction force we get
%
\begin{equation}
\begin{split}
\Fvec_i^h(t+\frac{h}{2}) &= - \frac{\gamma}{\alpha + \beta} \left( \vvec_i(t) - \uvec(\rvec_i,t) \right) - \frac{\alpha}{\alpha + \beta} \Fvec_i^\text{int} .
\end{split}
\end{equation}
%
%\begin{equation}\label{eq:modified-ebf}
%\begin{split}
%\Fvec_i^h(t+\frac{h}{2}) &= - \frac{m_i \rho a^3}{h(m_i + \rho a^3)} ( \vvec_i(t) - \uvec(\rvec_i,t) ) - \frac{\rho a^3}{m_i + \rho a^3} \Fvec_i^\text{int} \\
%&= - \frac{\rho a^3}{h(1 + \frac{\rho a^3}{m_i})} ( \vvec_i(t) - \uvec(\rvec_i,t) ) - \frac{1}{1 + \frac{m_i}{\rho a^3}} \Fvec_i^\text{int} \\
%&= - \frac{\gamma}{1+\alpha} ( \vvec_i(t) - \uvec(\rvec_i,t) ) -\frac{\alpha}{1+\alpha} \Fvec_i^\text{int} ,
%\end{split}
%\end{equation}
%
%The total force on the particle is hence
%%
%\begin{equation}
%\begin{split}
%\Fvec_i(t+\frac{h}{2}) &= \Fvec_i^h(t+\frac{h}{2}) + \Fvec_i^\text{int} \\ &= -
%\frac{\gamma}{\alpha + \beta} ( \vvec_i(t) - \uvec(\rvec_i,t) ) +
%\frac{\beta}{\alpha + \beta} \Fvec_i^\text{int} .
%\end{split}
%\end{equation}
%
The updates for the particle and flow velocity are thus
%
\begin{subequations}\label{eq:modified-ebf}
\begin{align}
\vvec_i(t+h) &= \vvec_i(t) 
- \frac{\alpha}{\alpha+\beta} \left( \vvec_i(t) - \uvec(\rvec_i,t) \right) 
+ \frac{h}{m_i} \frac{\beta}{\alpha+\beta} \Fvec_i^\text{int} , \\
\uvec(\rvec_i,t) &= \uvec(\rvec_i,t)
+ \frac{\beta}{\alpha+\beta} \left( \vvec_i(t) - \uvec(\rvec_i,t) \right) 
+ \frac{h}{\rho a^3} \frac{\alpha}{\alpha+\beta} \Fvec_i^\text{int} .
\end{align}
\end{subequations}
%
The modified EBF in Eq.~\eqref{eq:modified-ebf} leads to an
instantaneous no-slip boundary condition at the fluid-object
interface.  The implementation is straightforward and has the same
structure as the force coupling in Algorithm \ref{alg:fc}.
%
%This correction is necessary to ensure the no-slip boundary
%condition. 

If we introduce
%
\begin{align}\label{eq:immersion}
r &= \frac{\alpha}{\alpha + \beta} = \frac{1}{1 + \frac{m_i}{\rho
a^3}}, & (1-r) &= \frac{\beta}{\alpha+\beta} = \frac{1}{1 + \frac{\rho
a^3}{m_i}} ,
\end{align}
%
Equations~\eqref{eq:modified-ebf} can be interpreted as follows: A
fraction $r$ of the intrinsic force is included in the hydrodynamic
coupling and thus applied to the fluid, while the remaining fraction
$1-r$ is applied to the particle system. We remark that the same
equations are obtained, if we start out by assuming that the intrinsic
force is split between the fluid and the particles. That is,
Eqs.~\eqref{eq:modified-ebf} are a second order accurate scheme for
the differential equations
%
\begin{subequations}
\begin{align}
\frac{\d}{\d t} \vvec_i(t) &= - \frac{1}{m_i} \left[ \gamma \left( \vvec_i - \uvec(\rvec_i,t) \right) - \zetavec_i - (1-r) \Fvec_i^\text{int} \right] , \\
\frac{\d}{\d t} \uvec(\rvec_i,t) &= \frac{1}{\rho a^3} \left[ \gamma \left( \vvec_i - \uvec(\rvec_i,t) \right) - \zetavec_i  + r \Fvec_i^\text{int} \right] .
\end{align}
\end{subequations}
%
Equations~\eqref{eq:gamma-cond} and \eqref{eq:immersion} are then the
conditions under which the no-slip boundary condition is satisfied.

It is interesting to note that the distribution of the intrinsic force
is controlled by the ratio between the particle mass $m_i$ and
  the fluid mass within a unit cell $\rho a^3$, where $a$ is the
  lattice spacing. In the limit $m_i \gg \rho a^3$ we have $r
\rightarrow 0$, and the original EBF is recovered, where the intrinsic
force is applied to the particles only. Conversely, in the limit $m_i
\ll \rho a^3$ we have $r \rightarrow 1$ and recover the IBM where the
intrinsic force is applied to the fluid and the particles are simply
advected with the flow velocity. This shows that the immersed boundary
method corresponds to the limit where the particle mass becomes
negligible compared to the fluid mass. The dimensionless parameter $r$
controls whether the fluid is constrained to the particle velocity, or
whether the particles are constrained to the fluid flow.  We therefore
refer to $r$ as the ``immersion number''. It measures the relative
importance of inertia of the particles and the fluid.
%
The immersion number also depends on resolution: For fixed $\rho$ and
$m_i$, $r$ goes to zero in the continuum limit $a \rightarrow 0$.
%
% the i Moreover, it is dependent on the
%resolution: For fixed $\rho$ and $m_i$, the immersion number
%approaches one with increased grid spacing $a$. The surface points in
%turn represent a smaller portion of the mass of the immersed object
%and are more strongly constrained to the fluid flow.
%
For fixed resolution, $r$ depends only on the fraction $\alpha / \beta
= \rho a^3 / m_i$ and is independent of~$\gamma$. Hence $r$ can be
tuned to control the strength of the coupling with respect to the
transfer of intrinsic forces to the fluid. It is often desirable to
have a neutrally buoyant object, which is achieved when the immersion
number is $r=1/2$ or equivalently ${m_i}={\rho a^3}$. This choice is
also typical in other mesoscopic methods, for example, in coupling
elastic cells to multi-particle collision
dynamics~\cite{Peltomaki2013}.
%It is equivalent to $r=1/2$ which means that the intrinsic forces are
%split equally between the fluid and the particles.

The payoff for satisfying the no-slip boundary condition
instantaneously is the loss of freedom in tuning $\gamma$. The
physical friction and effective size of individual point-like
particles can not be tuned as easily any more.
%%
%\begin{equation}
%\begin{split}
%\frac{1}{\gamma_\text{eff}} &= \frac{1}{\gamma} + \frac{1}{6 \pi \eta g a}
%%= \frac{h}{2 \rho a^3 (1-r)} + \frac{1}{6 \pi \eta g a}
%= \frac{h}{m_i} + \frac{1}{6 \pi \eta g a} ,
%\end{split}
%\end{equation}
%%
%\begin{equation}
%\begin{split}
%\frac{1}{R_\text{eff}} &= \frac{6\pi\eta}{\gamma} + \frac{1}{ga} 
%%= \frac{6\pi\eta h}{2 \rho a^3 (1-r)} + \frac{1}{ga}
%= \frac{6 \pi \eta h}{m_i} + \frac{1}{ga} .
%\end{split} 
%\end{equation}
%%
However, this is fully compatible with the representation of the
solute object by a set of surface points: The object representation
has an explicit size, where the Lagrangian points are rather markers
than bead-like particles. The effective physical parameters of the
object are thus primarily determined by the boundary conditions on the
surface.



\section{Multiple time-step integration of force coupled systems}
\label{sec:multi-timestep}

Since the intrinsic forces are typically short-range, they often
require a time-step $\Delta t$ much smaller than a typical LB step
$h$. Reducing the overall time step can be costly, and as the dominant
computational effort is the fluid update, it is desirable to maintain
a large fluid time-step $h$. For accurate integration of the intrinsic
forces, we can sub-divide $h$ into an even integer number $2n$ of
sub-steps of size $\Delta t=h/(2n)$. Such a multiple time step scheme
can be written in terms of operator splitting as follows
\cite{Tuckerman1992}
%
\begin{equation}\label{eq:multistep}
\begin{split}
e^{h(\mathcal{S}+\mathcal{C}+\mathcal{P}+\mathcal{F})}
%&= \left[e^{\frac{h}{2}(\mathcal{P}+\mathcal{F})}\right] e^{h(\mathcal{S}+\mathcal{C})} \left[e^{\frac{h}{2}(\mathcal{P}+\mathcal{F})}\right] \\
&\approx \left[e^{n \Delta t (\mathcal{P}+\mathcal{F})}\right] e^{h(\mathcal{S}+\mathcal{C})} \left[e^{{n\Delta t}(\mathcal{P}+\mathcal{F})}\right] \\
&= \left[e^{\Delta t(\mathcal{P}+\mathcal{F})}\right]^{n} e^{h(\mathcal{S}+\mathcal{C})} \left[e^{\Delta t(\mathcal{P}+\mathcal{F})}\right]^{n} \\
&\approx \left[ \mathsf{P}^\frac{1}{4n} \mathsf{F}^\frac{1}{2n} \mathsf{P}^\frac{1}{4n} \right]^n [\mathsf{C}^\frac{1}{2}\mathsf{S}\mathsf{C}^\frac{1}{2}] \left[ \mathsf{P}^\frac{1}{4n} \mathsf{F}^\frac{1}{2n} \mathsf{P}^\frac{1}{4n} \right]^{n} . \\
\end{split}
\end{equation}
%
Equation~\eqref{eq:multistep} are two Verlet-like updates, each
iterated $n$ times with a time step $\Delta t$, and arranged
symmetrically around the LB update. The symmetric arrangement is
essential for accurate operator splitting, but it seems to receive
little attention in common software implementations.

Finally, we note that it is also possible to split the update of the
full system in the following way
\begin{equation}
e^{h(\mathcal{S}+\mathcal{C}+\mathcal{P}+\mathcal{F})}
\approx  \mathsf{F}_h^\frac12 \left[ \mathsf{P}^\frac{1}{4n} \mathsf{F}_\text{int}^\frac{1}{2n} \mathsf{P}^\frac{1}{4n} \right]^n [\mathsf{C}^\frac12\mathsf{S}\mathsf{C}^\frac12] \left[ \mathsf{P}^\frac{1}{4n} \mathsf{F}_\text{int}^\frac{1}{2n} \mathsf{P}^\frac{1}{4n} \right]^n \mathsf{F}_h^\frac12,
%e^{\frac{h}{2}(\mathcal{P}+\mathcal{F}_c)}[\mathsf{F}_h^\frac12 \mathsf{C}^\frac12\mathsf{S}\mathsf{C}^\frac12 \mathsf{F}_h^\frac12] e^{\frac{h}{2}(\mathcal{P}+\mathcal{F}_c)}
\end{equation}
where we have decomposed the force operator into the intrinsic part
$\mathsf{F}_\text{int}$ and the coupling part $\mathsf{F}_h$. In this
arrangement, the particle system can be integrated using an ODE solver
such as \textsf{LSODE} or \textsf{VODE}
\cite{Hindmarsh1983,Brown1989}. The fluid-particle interactions are
arranged at the outside in order to satisfy the no-slip boundary
condition after a full cycle through all parts. Moreover, the number
of evaluations of the hydrodynamic coupling force is reduced by a
factor $n$. The payoff is that this may lead to inaccuracies in
tracing the flow field by the particle positions.
%
A detailed numerical analysis of the accuracy, efficiency and
stability of the above splitting schemes is left for future work.



\section{Conclusions}

The operator splitting approach for the lattice Boltzmann method makes
it straightforward to derive accurate time-marching schemes for
various coupling methods. The proposed second-order schemes for the
immersed boundary method, force coupling and external boundary force
improve the accuracy compared to commonly implemented first-order
algorithms. The modified external boundary force unveils the close
relation of the immersed boundary method and force coupling.  All
three methods can be cast into the form of force coupling and thus
have a unified foundation. The coupling strength can be quantified by
the ratio between fluid and particle mass. In practice, this allows us
to control the coupling by appropriate choice of the immersion number.

While finite-time accuracy of the solutions may not be the main
concern if statistical properties in equilibrium or non-equilibrium
steady states are evaluated, accurate integration methods are still
favorable to maintain long-time stability
\cite{Shardlow2003,DeFabritiis2006,Bou-Rabee2010}. Splitting methods
are fairly established in molecular dynamics simulations and lead to
superior algorithms, in particular when applied to systems with
multiple scales and long-range forces
\cite{Tuckerman1992,McLachlan2002,Bou-Rabee2013}.
%
It is therefore anticipated that the fluid-particle coupling methods
discussed in this work will also benefit from a systematic operator
splitting approach.
%
%The improved coupling algorithms are expected to be useful
%in a wide range of applications, and they should become standard
%implementations in common software packages for the lattice Boltzmann
%method.



\section*{Acknowledgments}

I would like to thank Marisol Ripoll, Roland G. Winkler, Burkhard
D\"unweg and Anthony J. C. Ladd for valuable discussions and helpful
comments.
%
I also appreciate the careful and constructive review by an anonymous
referee.
%
Financial support from Volkswagen Stiftung under Grant No.~I/83576 is
gratefully acknowledged.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\appendix

\section{Spatial interpolation functions}
\label{sec:delta}

In the following, we outline the construction of the spatial
interpolation functions which are used in the hybrid coupling schemes
to map between the Eulerian grid and the Lagrangian particle system.
More details can be found, e.g., in \cite{Peskin2002}.

The three-dimensional delta function $\delta_a(\xvec)$ is taken
as a product of one-dimensional functions
%
\begin{equation}
\delta_a(\xvec) = \phi(\frac{x}{a}) \phi(\frac{y}{a}) \phi(\frac{z}{a}) ,
\end{equation}
%
where $x,y,z$ are the components of the vector $\xvec$. By postulate,
the function $\phi(r)$ should be continuous and satisfy the following
conditions \cite{Peskin2002}
%
\begin{subequations}
\begin{align}
%\phi(r) &= 0 \quad |r| \ge 2 \\
\sum_j \phi(r-j) &= 1 \quad \forall r \label{eq:delta1} , \\
\sum_j j \phi(r-j) &= r  \quad \forall r  \label{eq:delta2} , \\
\sum_j \left( \phi(r-j) \right)^2 &= C \quad \forall r \label{eq:delta3}  ,
\end{align}
\end{subequations}
%
where the sum runs over all integer numbers $j$.

For reasons of computational efficiency, $\phi(r)$ should also have
bounded support. Conditions \eqref{eq:delta1} and \eqref{eq:delta2}
imply
%
\begin{subequations}
\begin{align}
\sum_\xvec \delta_a(\xvec - \Xvec) &= 1 , \\ %\quad \forall \Xvec \\
\sum_\xvec \xvec \, \delta_a(\xvec - \Xvec) &= \Xvec ,%\quad \forall \Xvec \\
\end{align}
\end{subequations}
%
where the sum over $\xvec$ runs over all nodes of the Eulerian grid,
i.e., $\xvec = ( k a, l a, m a)$ and $k$, $l$, $m$ run over all
integer numbers.
%
These conditions guarantee that the total force and torque on the
Eulerian grid are the same as on the Lagrangian grid
%
\begin{subequations}
\begin{align}
\sum_\xvec a^3 \Fvec(\xvec,t) 
&= \sum_\xvec \sum_i \Fvec_i(t) \delta_a(\xvec - \Xvec_i) 
= \sum_i \Fvec_i(t) , \\
\sum_\xvec \xvec \times a^3 \Fvec(\xvec,t)
&= \sum_\xvec \sum_i \xvec \times \Fvec_i(t) \delta_a(\xvec - \Xvec_i)
= \sum_i \Xvec_i \times \Fvec_i(t) .
\end{align}
\end{subequations}
%
Moreover, \eqref{eq:delta1} and \eqref{eq:delta2} imply that smooth
functions are interpolated with second-order spatial accuracy.
%
The third condition \eqref{eq:delta3} can also be written as a Schwarz
inequality \cite{Peskin2002}
%
\begin{equation}\label{eq:delta4}
\left| \sum_j \phi(r_1-j) \phi(r_2-j) \right| \le C ,
\end{equation}
%
where the equal sign holds for $r_1=r_2$. This ensures that the
overlap of two Lagrangian particles is strongest when their positions
coincide, and the strength is independent of position. Condition
\eqref{eq:delta3} is thus a substitute for Galilean invariance, which
is impossible to satisfy if $\phi(r)$ has bounded support,
cf. \cite{Peskin2002}.

The three conditions \eqref{eq:delta1}-\eqref{eq:delta2} can be used to
determine a three-point interpolation function $\phi_3(r)$
%
\begin{equation}
\phi_3(r) = 
\begin{cases}
\frac{1}{3} \left( 1 + \sqrt{1 - 3r^2} \right) & 0 \le |r| \le \frac12 \\
\frac{1}{6} \left( 5 - 3|r| - \sqrt{6 |r| - 2 - 3r^2} \right) & \frac12 \le |r| \le \frac32 \\
0 & \frac32 \le |r| .
\end{cases}
\end{equation}
%
This function has a continuous derivative $\phi_3'(r)$, in contrast to
the linear interpolation used, e.g. in \cite{Ahlrichs1999}. It
therefore has the advantage that the flow gradient $\nabla \uvec$
varies smoothly and the particles do not experience jumps when
crossing grid lines. A numerical analysis of the properties of
different interpolation functions can be found in \cite{Duenweg2008}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\clearpage

\begin{algorithm*}
%
\caption{Collide-Stream-Collide Lattice Boltzmann}\label{alg:csc}
%
\begin{algorithmic}
  \Ensure{Parallelization has communicated processor boundaries}
  \Procedure{LBUpdate}{}
  %\Comment{$\gamma_k = 1 - \lambda_k$}
  %\Statex
  %\Comment{$\varphi_k = \sqrt{\mu b_k (1 - \gamma_k^2)}$}
  %\Statex
    \ForAll{$(x,y,z)$}
      \State \Call{LBPush}{$x,y,z$}
      \If{$(x,y,z) \in Interior$} \Call{LBRead}{$x-1,y-1,z-1$} \EndIf
    \EndFor
    \State $t \leftarrow t+h$
  \EndProcedure
\end{algorithmic}
%
\hrulefill
%
\begin{algorithmic}
  \Require{Random numbers $r_k(x,y,z)$}
  \Procedure{LBPush}{$x,y,z$}
    \ForAll{$k$}
      \State $m_k(x,y,z) \leftarrow m_k^\eq(x,y,z) + \gamma_k m_k^\nq(x,y,z) + \sqrt{\rho}\,\varphi_k r_k(x,y,z)$ 
      \State $f_i(x+c_{ix},y+c_{iy},z+c_{iz}) \leftarrow \sum_k a^{c_i}/b_k \cdot e_{ki} m_k(x,y,z) + G_i$ %\Comment{$G_i$: External force}
    \EndFor
  \EndProcedure
\end{algorithmic}
%
\hrulefill
%
\begin{algorithmic}
  \Require{Random numbers $r_k(x,y,z)$}
  \Procedure{LBRead}{$x,y,z$}
    \ForAll{$k$}
      \State $m_k(x,y,z) \leftarrow \sum_i e_{ki} f_i(x,y,z)$
      \State $m_k(x,y,z) \leftarrow m_k^\eq(x,y,z) + \gamma_k m_k^\nq(x,y,z) + \sqrt{\rho}\,\varphi_k r_k(x,y,z) + \sum_i e_{ki} G_i$
    \EndFor
  \EndProcedure
\end{algorithmic}
%
\end{algorithm*}

\begin{algorithm*}
%
\caption{Immersed Boundary Method}\label{alg:ibm}
%
\begin{algorithmic}
  \Procedure{ImmersedBoundaryMethod}{}
    \ForAll{$i \in Particles$}
      \State $\tilde\Xvec_i \leftarrow \Xvec_i + h/2 \cdot \Iop_a[\Xvec_i] \uvec$
      %\State $\Vvec_i \leftarrow \Iop_a[\tilde\Xvec_i] \uvec$
    \EndFor
    \State $\Fvec_i \leftarrow$ \Call{CalculateForces}{$\tilde\Xvec_i$}
    \ForAll{$i \in Particles$}
      \State $\Xvec_i \leftarrow \Xvec_i + h/2 \cdot \Iop_a[\tilde\Xvec_i] \uvec$
      \State \Call{SpreadForce}{$\Iop_a^*[\tilde\Xvec_i] \, \Fvec_i$}
    \EndFor
    \State \Call{LBUpdate}{}
    \ForAll{$i \in Particles$}
      \State \Call{SpreadForce}{$\Iop_a^*[\tilde\Xvec_i] \, \Fvec_i$}
      \State $\Xvec_i \leftarrow \Xvec_i + h/2 \cdot ( \Vvec_i + \Iop_a[\tilde\Xvec_i] \uvec )$
    \EndFor
  \EndProcedure
\end{algorithmic}
%
\end{algorithm*}



\begin{algorithm*}
%
\caption{Force Coupling}\label{alg:fc}
%
\begin{algorithmic}
  \Require{Random forces $\zetavec_i$}
%  \Comment{For no-EBF set $\gamma$ as in \eqref{eq:gamma-cond} and $r$ as in \eqref{eq:immersion}}
  \Procedure{CalculateCoupling}{$\rvec_i,\vvec_i,\Fvec_i,\uvec$}
    \ForAll{$i \in Particles$}
      \State $\Fvec_i^h \leftarrow - 1/(1+\alpha/2+\beta/2) \cdot ( \gamma (\vvec_i - \Iop_a[\rvec_i] \uvec ) - \zetavec_i + \alpha/2 \cdot \Fvec_i)$
    \EndFor
  \EndProcedure
%
\\\hrulefill
%
  \Procedure{ForceCoupling}{}
    \ForAll{$i \in Particles$}
      \State $\rvec_i \leftarrow \rvec_i + \Delta t/2 \cdot \vvec_i$
    \EndFor
    \State $\Fvec_i \leftarrow$ \Call{CalculateForces}{$\rvec_i$}
    \ForAll{$i \in Particles$}
      \State $\Fvec_i^h \leftarrow$ \Call{CalculateCoupling}{$\rvec_i,\vvec_i,\Fvec_i,\uvec$}
      \State $\vvec_i \leftarrow \vvec_i + \Delta t/2 \cdot (\Fvec_i^h + \Fvec_i)$
      \State \Call{SpreadForce}{$-\Iop_a^*[\rvec_i] \, \Fvec_i^h$}
    \EndFor
    \State \Call{LBUpdate}{}
    \ForAll{$i \in Particles$}
      \State $\Fvec_i^h \leftarrow$ \Call{CalculateCoupling}{$\rvec_i,\vvec_i,\Fvec_i,\uvec$}
      \State \Call{SpreadForce}{$-\Iop_a^*[\rvec_i] \, \Fvec_i^h$}
      \State $\vvec_i \leftarrow \vvec_i + \Delta t/2 \cdot (\Fvec_i^h + \Fvec_i)$
      \State $\rvec_i \leftarrow \rvec_i + \Delta t/2 \cdot \vvec_i$
    \EndFor       
  \EndProcedure
\end{algorithmic}
%
\end{algorithm*}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\clearpage

%\bibliographystyle{elsarticle-num}
\bibliographystyle{model1-num-names}
\bibliography{lb-operator-splitting}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}

lb-operator-splitting.bbl
\begin{thebibliography}{76}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\href}[2]{#2}
\providecommand{\path}[1]{#1}
\providecommand{\DOIprefix}{doi:}
\providecommand{\ArXivprefix}{arXiv:}
\providecommand{\URLprefix}{URL: }
\providecommand{\Pubmedprefix}{pmid:}
\providecommand{\doi}[1]{\href{http://dx.doi.org/#1}{\path{#1}}}
\providecommand{\Pubmed}[1]{\href{pmid:#1}{\path{#1}}}
\providecommand{\bibinfo}[2]{#2}
\ifx\xfnm\relax \def\xfnm[#1]{\unskip,\space#1}\fi
%Type = Book
\bibitem[{Succi(2001)}]{Succi-LB}
\bibinfo{author}{S.~Succi}, \bibinfo{title}{The Lattice {B}oltzmann Equation
  for Fluid Dynamics and Beyond}, \bibinfo{publisher}{Clarendon Press},
  \bibinfo{year}{2001}.
%Type = Article
\bibitem[{He and Luo(1997{\natexlab{a}})}]{He1997a}
\bibinfo{author}{X.~He}, \bibinfo{author}{L.-S. Luo},
\newblock \bibinfo{title}{A priori derivation of the lattice {B}oltzmann
  equation},
\newblock \bibinfo{journal}{Phys. Rev. E} \bibinfo{volume}{55}
  (\bibinfo{year}{1997}{\natexlab{a}}) \bibinfo{pages}{R6333--R6336}.
%Type = Article
\bibitem[{He and Luo(1997{\natexlab{b}})}]{He1997b}
\bibinfo{author}{X.~He}, \bibinfo{author}{L.-S. Luo},
\newblock \bibinfo{title}{Theory of the lattice {B}oltzmann method: {F}rom the
  {B}oltzmann equation to the lattice {B}oltzmann equation},
\newblock \bibinfo{journal}{Phys. Rev. E} \bibinfo{volume}{56}
  (\bibinfo{year}{1997}{\natexlab{b}}) \bibinfo{pages}{6811--6817}.
%Type = Article
\bibitem[{Shan and He(1998)}]{Shan1998}
\bibinfo{author}{X.~Shan}, \bibinfo{author}{X.~He},
\newblock \bibinfo{title}{Discretization of the velocity space in the solution
  of the {B}oltzmann equation},
\newblock \bibinfo{journal}{Phys. Rev. Lett.} \bibinfo{volume}{80}
  (\bibinfo{year}{1998}) \bibinfo{pages}{65--68}.
%Type = Article
\bibitem[{Junk et~al.(2005)Junk, Klar, and Luo}]{Junk2005a}
\bibinfo{author}{M.~Junk}, \bibinfo{author}{A.~Klar}, \bibinfo{author}{L.~S.
  Luo},
\newblock \bibinfo{title}{{A}symptotic analysis of the lattice {B}oltzmann
  equation},
\newblock \bibinfo{journal}{J. Comp. Phys.} \bibinfo{volume}{210}
  (\bibinfo{year}{2005}) \bibinfo{pages}{676--704}.
%Type = Article
\bibitem[{Krafczyk et~al.(2001)Krafczyk, T\"olke, Rank, and
  Schulz}]{Krafczyk2001}
\bibinfo{author}{M.~Krafczyk}, \bibinfo{author}{J.~T\"olke},
  \bibinfo{author}{E.~Rank}, \bibinfo{author}{M.~Schulz},
\newblock \bibinfo{title}{Two-dimensional simulation of fluid-structure
  interaction using lattice-{B}oltzmann methods},
\newblock \bibinfo{journal}{Computers \& Structures} \bibinfo{volume}{79}
  (\bibinfo{year}{2001}) \bibinfo{pages}{2031--2037}.
%Type = Incollection
\bibitem[{Geller et~al.(2006)Geller, T\"olke, and Krafczyk}]{Geller2006}
\bibinfo{author}{S.~Geller}, \bibinfo{author}{J.~T\"olke},
  \bibinfo{author}{M.~Krafczyk},
\newblock \bibinfo{title}{Lattice-{B}oltzmann method on quadtree-type grids for
  fluid-structure interaction},
\newblock in: \bibinfo{editor}{H.-J. Bungartz}, \bibinfo{editor}{M.~Sch\"afer}
  (Eds.), \bibinfo{booktitle}{Lecture Notes in Computational Science and
  Engineering}, volume~\bibinfo{volume}{53}, \bibinfo{publisher}{Springer},
  \bibinfo{year}{2006}, pp. \bibinfo{pages}{270--293}.
%Type = Article
\bibitem[{De~Rosis et~al.(2013)De~Rosis, Falcucci, Ubertini, and
  Ubertini}]{DeRosis2013}
\bibinfo{author}{A.~De~Rosis}, \bibinfo{author}{G.~Falcucci},
  \bibinfo{author}{S.~Ubertini}, \bibinfo{author}{F.~Ubertini},
\newblock \bibinfo{title}{A coupled lattice {B}oltzmann-finite element approach
  for two-dimensional fluid-structure interaction},
\newblock \bibinfo{journal}{Computers \& Fluids} \bibinfo{volume}{86}
  (\bibinfo{year}{2013}) \bibinfo{pages}{558--568}.
%Type = Article
\bibitem[{D\"unweg and Ladd(2008)}]{Duenweg2008}
\bibinfo{author}{B.~D\"unweg}, \bibinfo{author}{A.~J.~C. Ladd},
\newblock \bibinfo{title}{Lattice {B}oltzmann simulations of soft matter
  systems},
\newblock \bibinfo{journal}{Adv. Poly. Sci.} \bibinfo{volume}{221}
  (\bibinfo{year}{2008}) \bibinfo{pages}{89}.
%Type = Article
\bibitem[{Aidun and Clausen(2010)}]{Aidun2010}
\bibinfo{author}{C.~K. Aidun}, \bibinfo{author}{J.~R. Clausen},
\newblock \bibinfo{title}{Lattice-{B}oltzmann method for complex flows},
\newblock \bibinfo{journal}{Annu. Rev. Fluid Mech.} \bibinfo{volume}{42}
  (\bibinfo{year}{2010}) \bibinfo{pages}{439--472}.
%Type = Article
\bibitem[{Ladd(1994)}]{Ladd1994a}
\bibinfo{author}{A.~J.~C. Ladd},
\newblock \bibinfo{title}{Numerical simulations of particulate suspensions via
  a discretized {B}oltzmann equation. {P}art 1. {T}heoretical foundation},
\newblock \bibinfo{journal}{J. Fluid Mech.} \bibinfo{volume}{271}
  (\bibinfo{year}{1994}) \bibinfo{pages}{285--309}.
%Type = Article
\bibitem[{Ladd and Verberg(2001)}]{Ladd2001}
\bibinfo{author}{A.~J.~C. Ladd}, \bibinfo{author}{R.~Verberg},
\newblock \bibinfo{title}{Lattice-{B}oltzmann simulations of particle-fluid
  suspensions},
\newblock \bibinfo{journal}{J. Stat. Phys.} \bibinfo{volume}{104}
  (\bibinfo{year}{2001}) \bibinfo{pages}{1191--1251}.
%Type = Article
\bibitem[{Ahlrichs and D\"unweg(1999)}]{Ahlrichs1999}
\bibinfo{author}{P.~Ahlrichs}, \bibinfo{author}{B.~D\"unweg},
\newblock \bibinfo{title}{Simulation of a single polymer chain in solution by
  combining lattice {B}oltzmann and molecular dynamics},
\newblock \bibinfo{journal}{J. Chem. Phys.} \bibinfo{volume}{111}
  (\bibinfo{year}{1999}) \bibinfo{pages}{8225--8239}.
%Type = Article
\bibitem[{Pham et~al.(2009)Pham, Schiller, Prakash, and Duenweg}]{Pham2009}
\bibinfo{author}{T.~T. Pham}, \bibinfo{author}{U.~D. Schiller},
  \bibinfo{author}{J.~R. Prakash}, \bibinfo{author}{B.~Duenweg},
\newblock \bibinfo{title}{Implicit and explicit solvent models for the
  simulation of a single polymer chain in solution: {L}attice {B}oltzmann vs
  {B}rownian dynamics},
\newblock \bibinfo{journal}{J. Chem. Phys.} \bibinfo{volume}{131}
  (\bibinfo{year}{2009}) \bibinfo{pages}{164114}.
%Type = Article
\bibitem[{Ladd et~al.(2009)Ladd, Kekre, and Butler}]{Ladd2009}
\bibinfo{author}{A.~J.~C. Ladd}, \bibinfo{author}{R.~Kekre},
  \bibinfo{author}{J.~E. Butler},
\newblock \bibinfo{title}{Comparison of the static and dynamic properties of a
  semiflexible polymer using lattice {B}oltzmann and {B}rownian-dynamics
  simulations},
\newblock \bibinfo{journal}{Phys. Rev. E} \bibinfo{volume}{80}
  (\bibinfo{year}{2009}) \bibinfo{pages}{036704}.
%Type = Article
\bibitem[{Alexeev et~al.(2006)Alexeev, Verberg, and Balazs}]{Alexeev2005}
\bibinfo{author}{A.~Alexeev}, \bibinfo{author}{R.~Verberg},
  \bibinfo{author}{A.~C. Balazs},
\newblock \bibinfo{title}{{M}odeling the interactions between deformable
  capsules rolling on a compliant surface},
\newblock \bibinfo{journal}{Soft Matter} \bibinfo{volume}{2}
  (\bibinfo{year}{2006}) \bibinfo{pages}{499--509}.
%Type = Article
\bibitem[{Dupin et~al.(2007)Dupin, Halliday, Care, Lyuba, and Munn}]{Dupin2007}
\bibinfo{author}{M.~M. Dupin}, \bibinfo{author}{I.~Halliday},
  \bibinfo{author}{C.~M. Care}, \bibinfo{author}{A.~Lyuba},
  \bibinfo{author}{L.~L. Munn},
\newblock \bibinfo{title}{Modeling the flow of dense suspensions of deformable
  particles in three dimensions},
\newblock \bibinfo{journal}{Phys. Rev. E} \bibinfo{volume}{75}
  (\bibinfo{year}{2007}) \bibinfo{pages}{066707}.
%Type = Article
\bibitem[{Wu and Aidun(2010)}]{Wu2010}
\bibinfo{author}{J.~Wu}, \bibinfo{author}{C.~K. Aidun},
\newblock \bibinfo{title}{Simulating {3D} deformable particle suspensions using
  lattice {B}oltzmann method with discrete external boundary force},
\newblock \bibinfo{journal}{Int. J. Numer. Meth. Fluids} \bibinfo{volume}{62}
  (\bibinfo{year}{2010}) \bibinfo{pages}{765--783}.
%Type = Article
\bibitem[{Brownlee et~al.(2007)Brownlee, Gorban, and Levesley}]{Brownlee2007}
\bibinfo{author}{R.~A. Brownlee}, \bibinfo{author}{A.~N. Gorban},
  \bibinfo{author}{J.~Levesley},
\newblock \bibinfo{title}{Stability and stabilization of the lattice
  {B}oltzmann method},
\newblock \bibinfo{journal}{Phys. Rev. E} \bibinfo{volume}{75}
  (\bibinfo{year}{2007}) \bibinfo{pages}{036711}.
%Type = Article
\bibitem[{Peskin(2002)}]{Peskin2002}
\bibinfo{author}{C.~S. Peskin},
\newblock \bibinfo{title}{The immersed boundary method},
\newblock \bibinfo{journal}{Acta Numer.} \bibinfo{volume}{11}
  (\bibinfo{year}{2002}) \bibinfo{pages}{479--517}.
%Type = Article
\bibitem[{Serrano et~al.(2006)Serrano, De~Fabritiis, Espa\~{n}ol, and
  Coveney}]{Serrano2006}
\bibinfo{author}{M.~Serrano}, \bibinfo{author}{G.~De~Fabritiis},
  \bibinfo{author}{P.~Espa\~{n}ol}, \bibinfo{author}{P.~Coveney},
\newblock \bibinfo{title}{A stochastic {T}rotter integration scheme for
  dissipative particle dynamics},
\newblock \bibinfo{journal}{Math. Comp. Sim.} \bibinfo{volume}{72}
  (\bibinfo{year}{2006}) \bibinfo{pages}{190--194}.
%Type = Article
\bibitem[{Schlick et~al.(1998)Schlick, Mandziuk, Skeel, and
  Srinivas}]{Schlick1998}
\bibinfo{author}{T.~Schlick}, \bibinfo{author}{M.~Mandziuk},
  \bibinfo{author}{R.~D. Skeel}, \bibinfo{author}{K.~Srinivas},
\newblock \bibinfo{title}{Nonlinear resonance artifacts in molecular dynamics
  simulations},
\newblock \bibinfo{journal}{J. Comp. Phys.} \bibinfo{volume}{140}
  (\bibinfo{year}{1998}) \bibinfo{pages}{1--29}.
%Type = Article
\bibitem[{Tuckerman et~al.(1992)Tuckerman, Berne, and Martyna}]{Tuckerman1992}
\bibinfo{author}{M.~Tuckerman}, \bibinfo{author}{B.~J. Berne},
  \bibinfo{author}{G.~J. Martyna},
\newblock \bibinfo{title}{Reversible multiple time scale molecular dynamics},
\newblock \bibinfo{journal}{J. Chem. Phys.} \bibinfo{volume}{97}
  (\bibinfo{year}{1992}) \bibinfo{pages}{1990--2001}.
%Type = Article
\bibitem[{Ricci and Ciccotti(2003)}]{Ricci2003}
\bibinfo{author}{A.~Ricci}, \bibinfo{author}{G.~Ciccotti},
\newblock \bibinfo{title}{Algorithms for {B}rownian dynamics},
\newblock \bibinfo{journal}{Mol. Phys.} \bibinfo{volume}{101}
  (\bibinfo{year}{2003}) \bibinfo{pages}{1927--1931}.
%Type = Article
\bibitem[{Bussi and Parrinello(2007)}]{Bussi2007}
\bibinfo{author}{G.~Bussi}, \bibinfo{author}{M.~Parrinello},
\newblock \bibinfo{title}{Accurate sampling using {L}angevin dynamics},
\newblock \bibinfo{journal}{Phys. Rev. E} \bibinfo{volume}{75}
  (\bibinfo{year}{2007}) \bibinfo{pages}{056707}.
%Type = Article
\bibitem[{Trotter(1959)}]{Trotter1959}
\bibinfo{author}{H.~F. Trotter},
\newblock \bibinfo{title}{On the product of semi-groups of operators},
\newblock \bibinfo{journal}{Proc. Amer. Math. Soc.} \bibinfo{volume}{10}
  (\bibinfo{year}{1959}) \bibinfo{pages}{545--551}.
%Type = Incollection
\bibitem[{Kato(1978)}]{Kato1978}
\bibinfo{author}{T.~Kato},
\newblock \bibinfo{title}{Trotter's product formula for an arbitrary pair of
  self-adjoint contraction semigroups},
\newblock in: \bibinfo{booktitle}{Topics in functional analysis (essays
  dedicated to {M}. {G}. {K}re\u\i n on the occasion of his 70th birthday)},
  volume~\bibinfo{volume}{3} of \textit{\bibinfo{series}{Adv. in Math. Suppl.
  Stud.}}, \bibinfo{publisher}{Academic Press}, \bibinfo{address}{New York},
  \bibinfo{year}{1978}, pp. \bibinfo{pages}{185--195}.
%Type = Article
\bibitem[{Suzuki(1985)}]{Suzuki1985}
\bibinfo{author}{M.~Suzuki},
\newblock \bibinfo{title}{Decomposition formulas of exponential operators and
  {L}ie exponentials with some applications to quantum mechanics and
  statistical physics},
\newblock \bibinfo{journal}{J. Math. Phys.} \bibinfo{volume}{26}
  (\bibinfo{year}{1985}) \bibinfo{pages}{601--612}.
%Type = Article
\bibitem[{McLachlan and Quispel(2002)}]{McLachlan2002}
\bibinfo{author}{R.~I. McLachlan}, \bibinfo{author}{G.~R.~W. Quispel},
\newblock \bibinfo{title}{Splitting methods},
\newblock \bibinfo{journal}{Acta Numer.} \bibinfo{volume}{11}
  (\bibinfo{year}{2002}) \bibinfo{pages}{341--434}.
%Type = Article
\bibitem[{Bou-Rabee(2013)}]{Bou-Rabee2013}
\bibinfo{author}{N.~Bou-Rabee},
\newblock \bibinfo{title}{Time integrators for molecular dynamics},
\newblock \bibinfo{journal}{Entropy} \bibinfo{volume}{16}
  (\bibinfo{year}{2013}) \bibinfo{pages}{138--162}.
%Type = Article
\bibitem[{Dellar(2013)}]{Dellar2013}
\bibinfo{author}{P.~J. Dellar},
\newblock \bibinfo{title}{An interpretation and derivation of the lattice
  {B}oltzmann method using {S}trang splitting},
\newblock \bibinfo{journal}{Comput. Math. Appl.} \bibinfo{volume}{65}
  (\bibinfo{year}{2013}) \bibinfo{pages}{129--141}.
%Type = Book
\bibitem[{Frenkel and Smit(2001)}]{FrenkelSmit}
\bibinfo{author}{D.~Frenkel}, \bibinfo{author}{B.~Smit},
  \bibinfo{title}{Understanding Molecular Simulation: From Algorithms to
  Applications}, \bibinfo{publisher}{Academic Press}, \bibinfo{year}{2001}.
%Type = Article
\bibitem[{Skeel et~al.(1997)Skeel, Zhang, and Schlick}]{Skeel1997}
\bibinfo{author}{R.~Skeel}, \bibinfo{author}{G.~Zhang},
  \bibinfo{author}{T.~Schlick},
\newblock \bibinfo{title}{A family of symplectic integrators: Stability,
  accuracy, and molecular dynamics applications},
\newblock \bibinfo{journal}{SIAM J. Sci. Comput.} \bibinfo{volume}{18}
  (\bibinfo{year}{1997}) \bibinfo{pages}{203--222}.
%Type = Article
\bibitem[{Shardlow(2003)}]{Shardlow2003}
\bibinfo{author}{T.~Shardlow},
\newblock \bibinfo{title}{Splitting for dissipative particle dynamics},
\newblock \bibinfo{journal}{SIAM J. Sci. Comput.} \bibinfo{volume}{24}
  (\bibinfo{year}{2003}) \bibinfo{pages}{1267--1282}.
%Type = Article
\bibitem[{De~Fabritiis et~al.(2006)De~Fabritiis, Serrano, Espa\~nol, and
  Coveney}]{DeFabritiis2006}
\bibinfo{author}{G.~De~Fabritiis}, \bibinfo{author}{M.~Serrano},
  \bibinfo{author}{P.~Espa\~nol}, \bibinfo{author}{P.~Coveney},
\newblock \bibinfo{title}{Efficient numerical integrators for stochastic
  models},
\newblock \bibinfo{journal}{Physica A} \bibinfo{volume}{361}
  (\bibinfo{year}{2006}) \bibinfo{pages}{429--440}.
%Type = Article
\bibitem[{Malevanets and Kapral(1999)}]{Malevanets1999}
\bibinfo{author}{A.~Malevanets}, \bibinfo{author}{R.~Kapral},
\newblock \bibinfo{title}{Mesoscopic model for solvent dynamics},
\newblock \bibinfo{journal}{J. Chem. Phys.} \bibinfo{volume}{110}
  (\bibinfo{year}{1999}) \bibinfo{pages}{8605--8613}.
%Type = Article
\bibitem[{Gompper et~al.(2009)Gompper, Ihle, Kroll, and Winkler}]{Gompper2009}
\bibinfo{author}{G.~Gompper}, \bibinfo{author}{T.~Ihle}, \bibinfo{author}{D.~M.
  Kroll}, \bibinfo{author}{R.~G. Winkler},
\newblock \bibinfo{title}{Multi-particle collision dynamics: A particle-based
  mesoscale simulation approach to the hydrodynamics of complex fluids},
\newblock \bibinfo{journal}{Adv. Poly. Sci.} \bibinfo{volume}{221}
  (\bibinfo{year}{2009}) \bibinfo{pages}{1--87}.
%Type = Book
\bibitem[{Cercignani(1988)}]{Cercignani1988}
\bibinfo{author}{C.~Cercignani}, \bibinfo{title}{The {B}oltzmann Equation and
  Its Applications}, volume~\bibinfo{volume}{67} of
  \textit{\bibinfo{series}{Applied Mathematical Sciences}},
  \bibinfo{publisher}{Springer}, \bibinfo{year}{1988}.
%Type = Article
\bibitem[{Bhatnagar et~al.(1954)Bhatnagar, Gross, and Krook}]{Bhatnagar1954}
\bibinfo{author}{P.~L. Bhatnagar}, \bibinfo{author}{E.~P. Gross},
  \bibinfo{author}{M.~Krook},
\newblock \bibinfo{title}{A model for collision processes in gases. i. small
  amplitude processes in charged and neutral one-component systems},
\newblock \bibinfo{journal}{Phys. Rev.} \bibinfo{volume}{94}
  (\bibinfo{year}{1954}) \bibinfo{pages}{511--525}.
%Type = Article
\bibitem[{Gross and Jackson(1959)}]{Gross1959}
\bibinfo{author}{E.~P. Gross}, \bibinfo{author}{E.~A. Jackson},
\newblock \bibinfo{title}{Kinetic models and the linearized {B}oltzmann
  equation},
\newblock \bibinfo{journal}{Phys. Fluids} \bibinfo{volume}{2}
  (\bibinfo{year}{1959}) \bibinfo{pages}{432--441}.
%Type = Article
\bibitem[{Shan et~al.(2006)Shan, Yuan, and Chen}]{Shan2006}
\bibinfo{author}{X.~Shan}, \bibinfo{author}{X.-F. Yuan},
  \bibinfo{author}{H.~Chen},
\newblock \bibinfo{title}{Kinetic theory representation of hydrodynamics: a way
  beyond the {N}avier-{S}tokes equation},
\newblock \bibinfo{journal}{J. Fluid Mech.} \bibinfo{volume}{550}
  (\bibinfo{year}{2006}) \bibinfo{pages}{413--441}.
%Type = Article
\bibitem[{Broadwell(1964)}]{Broadwell1964a}
\bibinfo{author}{J.~E. Broadwell},
\newblock \bibinfo{title}{Study of rarefied shear flow by the discrete velocity
  method},
\newblock \bibinfo{journal}{J. Fluid Mech.} \bibinfo{volume}{19}
  (\bibinfo{year}{1964}) \bibinfo{pages}{401--414}.
%Type = Article
\bibitem[{Luo(2000)}]{Luo2000}
\bibinfo{author}{L.-S. Luo},
\newblock \bibinfo{title}{Some recent results on discrete velocity models and
  ramifications for lattice {B}oltzmann equation},
\newblock \bibinfo{journal}{Comp. Phys. Comm.} \bibinfo{volume}{129}
  (\bibinfo{year}{2000}) \bibinfo{pages}{63--74}.
%Type = Article
\bibitem[{He et~al.(1998)He, Shan, and Doolen}]{He1998}
\bibinfo{author}{X.~He}, \bibinfo{author}{X.~Shan}, \bibinfo{author}{G.~D.
  Doolen},
\newblock \bibinfo{title}{Discrete {B}oltzmann equation model for nonideal
  gases},
\newblock \bibinfo{journal}{Phys. Rev. E} \bibinfo{volume}{57}
  (\bibinfo{year}{1998}) \bibinfo{pages}{R13--R16}.
%Type = Article
\bibitem[{Benzi et~al.(1992)Benzi, Succi, and Vergalossa}]{Benzi1992}
\bibinfo{author}{R.~Benzi}, \bibinfo{author}{S.~Succi},
  \bibinfo{author}{M.~Vergalossa},
\newblock \bibinfo{title}{The lattice {B}oltzmann equation: theory and
  applications},
\newblock \bibinfo{journal}{Phys. Rep.} \bibinfo{volume}{222}
  (\bibinfo{year}{1992}) \bibinfo{pages}{145--197}.
%Type = Article
\bibitem[{Chen and Doolen(1998)}]{Chen1998}
\bibinfo{author}{S.~Chen}, \bibinfo{author}{G.~D. Doolen},
\newblock \bibinfo{title}{Lattice {B}oltzmann method for fluid flows},
\newblock \bibinfo{journal}{Annu. Rev. Fluid Mech.} \bibinfo{volume}{30}
  (\bibinfo{year}{1998}) \bibinfo{pages}{329--364}.
%Type = Article
\bibitem[{He et~al.(1998)He, Chen, and Doolen}]{He1998a}
\bibinfo{author}{X.~He}, \bibinfo{author}{S.~Chen}, \bibinfo{author}{G.~D.
  Doolen},
\newblock \bibinfo{title}{A novel thermal model for the lattice {B}oltzmann
  method in incompressible limit},
\newblock \bibinfo{journal}{J. Comp. Phys.} \bibinfo{volume}{146}
  (\bibinfo{year}{1998}) \bibinfo{pages}{282--300}.
%Type = Article
\bibitem[{Guo et~al.(2002)Guo, Zheng, and Shi}]{Guo2002}
\bibinfo{author}{Z.~Guo}, \bibinfo{author}{C.~Zheng}, \bibinfo{author}{B.~Shi},
\newblock \bibinfo{title}{Discrete lattice effects on the forcing term in the
  lattice {B}oltzmann method},
\newblock \bibinfo{journal}{Phys. Rev. E} \bibinfo{volume}{65}
  (\bibinfo{year}{2002}) \bibinfo{pages}{046308}.
%Type = Article
\bibitem[{Dellar(2003)}]{Dellar2003}
\bibinfo{author}{P.~J. Dellar},
\newblock \bibinfo{title}{Incompressible limits of lattice {B}oltzmann
  equations using multiple relaxation times},
\newblock \bibinfo{journal}{J. Comp. Phys.} \bibinfo{volume}{190}
  (\bibinfo{year}{2003}) \bibinfo{pages}{351--370}.
%Type = Article
\bibitem[{H\'enon(1987)}]{Henon1987}
\bibinfo{author}{M.~H\'enon},
\newblock \bibinfo{title}{Viscosity of a lattice gas},
\newblock \bibinfo{journal}{Complex Sys.} \bibinfo{volume}{1}
  (\bibinfo{year}{1987}) \bibinfo{pages}{763--789}.
%Type = Article
\bibitem[{Xu and He(2003)}]{Xu2003}
\bibinfo{author}{K.~Xu}, \bibinfo{author}{X.~He},
\newblock \bibinfo{title}{Lattice {B}oltzmann method and gas-kinetic {BGK}
  scheme in the low-{M}ach number viscous flow simulations},
\newblock \bibinfo{journal}{J. Comp. Phys.} \bibinfo{volume}{190}
  (\bibinfo{year}{2003}) \bibinfo{pages}{100--117}.
%Type = Article
\bibitem[{Ohwada(1998)}]{Ohwada1998}
\bibinfo{author}{T.~Ohwada},
\newblock \bibinfo{title}{Higher order approximation methods for the
  {B}oltzmann equation},
\newblock \bibinfo{journal}{J. Comp. Phys.} \bibinfo{volume}{139}
  (\bibinfo{year}{1998}) \bibinfo{pages}{1--14}.
%Type = Article
\bibitem[{Strang(1968)}]{Strang1968}
\bibinfo{author}{G.~Strang},
\newblock \bibinfo{title}{On the construction and comparison of difference
  schemes},
\newblock \bibinfo{journal}{SIAM J. Numer. Anal.} \bibinfo{volume}{5}
  (\bibinfo{year}{1968}) \bibinfo{pages}{506--517}.
%Type = Book
\bibitem[{Gr\"obner(1960)}]{Groebner1960}
\bibinfo{author}{W.~Gr\"obner}, \bibinfo{title}{Die {L}ie-Reihen und ihre
  {A}nwendungen}, \bibinfo{publisher}{Deutscher Verlag der Wissenschaften},
  \bibinfo{year}{1960}.
%Type = Book
\bibitem[{Varadarajan(1974)}]{Varadarajan1974}
\bibinfo{author}{V.~S. Varadarajan}, \bibinfo{title}{Lie groups, Lie algebras,
  and their representations}, \bibinfo{publisher}{Prentice-Hall},
  \bibinfo{year}{1974}.
%Type = Inbook
\bibitem[{Hairer et~al.(2006)Hairer, Lubich, and Wanner}]{Hairer2006}
\bibinfo{author}{E.~Hairer}, \bibinfo{author}{C.~Lubich},
  \bibinfo{author}{G.~Wanner}, \bibinfo{title}{Geometric Numerical
  Integration}, \bibinfo{edition}{2nd} ed., \bibinfo{publisher}{Springer},
  \bibinfo{year}{2006}, pp. \bibinfo{pages}{51--96}.
%Type = Article
\bibitem[{Yoshida(1990)}]{Yoshida1990}
\bibinfo{author}{H.~Yoshida},
\newblock \bibinfo{title}{Construction of higher order symplectic integrators},
\newblock \bibinfo{journal}{Physics Letters A} \bibinfo{volume}{150}
  (\bibinfo{year}{1990}) \bibinfo{pages}{262 -- 268}.
%Type = Article
\bibitem[{Bobylev and Ohwada(2001)}]{Bobylev2001}
\bibinfo{author}{A.~Bobylev}, \bibinfo{author}{T.~Ohwada},
\newblock \bibinfo{title}{The error of the splitting scheme for solving
  evolutionary equations},
\newblock \bibinfo{journal}{Appl. Math. Lett.} \bibinfo{volume}{14}
  (\bibinfo{year}{2001}) \bibinfo{pages}{45--48}.
%Type = Article
\bibitem[{Salmon(1999{\natexlab{a}})}]{Salmon1999a}
\bibinfo{author}{R.~Salmon},
\newblock \bibinfo{title}{The lattice {B}oltzmann method as a basis for ocean
  circulation modeling},
\newblock \bibinfo{journal}{J. Marine Res.} \bibinfo{volume}{57}
  (\bibinfo{year}{1999}{\natexlab{a}}) \bibinfo{pages}{503--535}.
%Type = Article
\bibitem[{Salmon(1999{\natexlab{b}})}]{Salmon1999b}
\bibinfo{author}{R.~Salmon},
\newblock \bibinfo{title}{Lattice {B}oltzmann solutions of the
  three-dimensional planetary geostrophic equations},
\newblock \bibinfo{journal}{J. Marine Res.} \bibinfo{volume}{57}
  (\bibinfo{year}{1999}{\natexlab{b}}) \bibinfo{pages}{847--884}.
%Type = Article
\bibitem[{Adhikari et~al.(2005)Adhikari, Stratford, Cates, and
  Wagner}]{Adhikari2005}
\bibinfo{author}{R.~Adhikari}, \bibinfo{author}{K.~Stratford},
  \bibinfo{author}{M.~E. Cates}, \bibinfo{author}{A.~J. Wagner},
\newblock \bibinfo{title}{Fluctuating lattice {B}oltzmann},
\newblock \bibinfo{journal}{Europhys. Lett.} \bibinfo{volume}{71}
  (\bibinfo{year}{2005}) \bibinfo{pages}{473}.
%Type = Article
\bibitem[{D\"unweg et~al.(2007)D\"unweg, Schiller, and Ladd}]{Duenweg2007}
\bibinfo{author}{B.~D\"unweg}, \bibinfo{author}{U.~D. Schiller},
  \bibinfo{author}{A.~J.~C. Ladd},
\newblock \bibinfo{title}{Statistical mechanics of the fluctuating lattice
  {B}oltzmann equation},
\newblock \bibinfo{journal}{Phys. Rev. E} \bibinfo{volume}{76}
  (\bibinfo{year}{2007}) \bibinfo{pages}{036704}.
%Type = Article
\bibitem[{Luo(1998)}]{Luo1998}
\bibinfo{author}{L.-S. Luo},
\newblock \bibinfo{title}{Unified theory of lattice {B}oltzmann models for
  nonideal gases},
\newblock \bibinfo{journal}{Phys. Rev. Lett.} \bibinfo{volume}{81}
  (\bibinfo{year}{1998}) \bibinfo{pages}{1618--1621}.
%Type = Article
\bibitem[{Martys et~al.(1998)Martys, Shan, and Chen}]{Martys1998}
\bibinfo{author}{N.~S. Martys}, \bibinfo{author}{X.~Shan},
  \bibinfo{author}{H.~Chen},
\newblock \bibinfo{title}{Evaluation of the external force term in the discrete
  {B}oltzmann equation},
\newblock \bibinfo{journal}{Phys. Rev. E} \bibinfo{volume}{58}
  (\bibinfo{year}{1998}) \bibinfo{pages}{6855--6857}.
%Type = Article
\bibitem[{Ubertini et~al.(2010)Ubertini, Asinari, and Succi}]{Ubertini2010}
\bibinfo{author}{S.~Ubertini}, \bibinfo{author}{P.~Asinari},
  \bibinfo{author}{S.~Succi},
\newblock \bibinfo{title}{Three ways to lattice {B}boltzmann: {A} unified
  time-marching picture},
\newblock \bibinfo{journal}{Phys. Rev. E} \bibinfo{volume}{81}
  (\bibinfo{year}{2010}) \bibinfo{pages}{016311}.
%Type = Book
\bibitem[{Chapman and Cowling(1990)}]{ChapmanCowling}
\bibinfo{author}{S.~Chapman}, \bibinfo{author}{T.~G. Cowling},
  \bibinfo{title}{The Mathematical Theory of Non-Uniform Gases}, Cambridge
  Mathematical Library Series, \bibinfo{edition}{3rd} ed.,
  \bibinfo{publisher}{Cambridge University Press}, \bibinfo{year}{1990}.
%Type = Article
\bibitem[{McNamara and Alder(1993)}]{McNamara1993}
\bibinfo{author}{G.~McNamara}, \bibinfo{author}{B.~Alder},
\newblock \bibinfo{title}{Analysis of the lattice {B}oltzmann treatment of
  hydrodynamics},
\newblock \bibinfo{journal}{Physica A} \bibinfo{volume}{194}
  (\bibinfo{year}{1993}) \bibinfo{pages}{218--228}.
%Type = Article
\bibitem[{Dellar(2001)}]{Dellar2001}
\bibinfo{author}{P.~J. Dellar},
\newblock \bibinfo{title}{Bulk and shear viscosities in lattice {B}oltzmann
  equations},
\newblock \bibinfo{journal}{Phys. Rev. E} \bibinfo{volume}{64}
  (\bibinfo{year}{2001}) \bibinfo{pages}{031203}.
%Type = Article
\bibitem[{Wittmann et~al.(2013)Wittmann, Zeiser, Hager, and
  Wellein}]{Wittmann2013}
\bibinfo{author}{M.~Wittmann}, \bibinfo{author}{T.~Zeiser},
  \bibinfo{author}{G.~Hager}, \bibinfo{author}{G.~Wellein},
\newblock \bibinfo{title}{Comparison of different propagation steps for lattice
  {B}oltzmann methods},
\newblock \bibinfo{journal}{Comput. Math. Appl.} \bibinfo{volume}{65}
  (\bibinfo{year}{2013}) \bibinfo{pages}{924--935}.
%Type = Article
\bibitem[{Lai and Peskin(2000)}]{Lai2000}
\bibinfo{author}{M.-C. Lai}, \bibinfo{author}{C.~S. Peskin},
\newblock \bibinfo{title}{An immersed boundary method with formal second-order
  accuracy and reduced numerical viscosity},
\newblock \bibinfo{journal}{J. Comp. Phys.} \bibinfo{volume}{160}
  (\bibinfo{year}{2000}) \bibinfo{pages}{705--719}.
%Type = Article
\bibitem[{Mori and Peskin(2008)}]{Mori2008}
\bibinfo{author}{Y.~Mori}, \bibinfo{author}{C.~S. Peskin},
\newblock \bibinfo{title}{Implicit second-order immersed boundary methods with
  boundary mass},
\newblock \bibinfo{journal}{Comput. Meth. Appl. Mech. Eng.}
  \bibinfo{volume}{197} (\bibinfo{year}{2008}) \bibinfo{pages}{2049--2067}.
%Type = Article
\bibitem[{Goldstein et~al.(1993)Goldstein, Handler, and
  Sirovich}]{Goldstein1993}
\bibinfo{author}{D.~Goldstein}, \bibinfo{author}{R.~Handler},
  \bibinfo{author}{L.~Sirovich},
\newblock \bibinfo{title}{Modeling a no-slip flow boundary with an external
  force field},
\newblock \bibinfo{journal}{J. Comp. Phys.} \bibinfo{volume}{105}
  (\bibinfo{year}{1993}) \bibinfo{pages}{354--366}.
%Type = Article
\bibitem[{Peltom\"aki and Gompper(2013)}]{Peltomaki2013}
\bibinfo{author}{M.~Peltom\"aki}, \bibinfo{author}{G.~Gompper},
\newblock \bibinfo{title}{Sedimentation of single red blood cells},
\newblock \bibinfo{journal}{Soft Matter} \bibinfo{volume}{9}
  (\bibinfo{year}{2013}) \bibinfo{pages}{8346--8358}.
%Type = Article
\bibitem[{Hindmarsh(1983)}]{Hindmarsh1983}
\bibinfo{author}{A.~C. Hindmarsh},
\newblock \bibinfo{title}{{ODEPACK}, a systematized collection of {ODE}
  solvers},
\newblock \bibinfo{journal}{IMACS Trans. Sci. Comput.} \bibinfo{volume}{1}
  (\bibinfo{year}{1983}) \bibinfo{pages}{55--64}.
%Type = Article
\bibitem[{Brown et~al.(1989)Brown, Byrne, and Hindmarsh}]{Brown1989}
\bibinfo{author}{P.~Brown}, \bibinfo{author}{G.~Byrne},
  \bibinfo{author}{A.~Hindmarsh},
\newblock \bibinfo{title}{{VODE}: A variable-coefficient {ODE} solver},
\newblock \bibinfo{journal}{SIAM J. Sci. Stat. Comput.} \bibinfo{volume}{10}
  (\bibinfo{year}{1989}) \bibinfo{pages}{1038--1051}.
%Type = Article
\bibitem[{Bou-Rabee and Owhadi(2010)}]{Bou-Rabee2010}
\bibinfo{author}{N.~Bou-Rabee}, \bibinfo{author}{H.~Owhadi},
\newblock \bibinfo{title}{Long-run accuracy of variational integrators in the
  stochastic context},
\newblock \bibinfo{journal}{SIAM J. Numer. Anal.} \bibinfo{volume}{48}
  (\bibinfo{year}{2010}) \bibinfo{pages}{278--297}.

\end{thebibliography}
