remarkRemark hypothesisHypothesis claimClaim MRT-LB model based four-level FD
scheme for 1D diffusion equationY.X. Lin, N. Hong, B.C. Shi and Z. H. Chai

# A multiple-relaxation-time lattice Boltzmann model based four-level finite-
difference scheme for one-dimensional diffusion equation††thanks: Submitted to
the editors DATE. This work was financially supported by the National Natural
Science Foundation of China (Grants No. 12072127 and No. 51836003).

Yuxin Lin School of Mathematics and Statistics, Huazhong University of Science
and Technology, Wuhan, 430074, China.  Ning Hong School of General Education,
Wuchang University of Technology, Wuhan 430223, China. Baochang Shi School of
Mathematics and Statistics, Huazhong University of Science and Technology,
Wuhan, 430074, China; Hubei Key Laboratory of Engineering Modeling and
Scientific Computing, Huazhong University of Science and Technology, Wuhan
430074, China. Zhenhua Chai22footnotemark: 2 Corresponding author. School of
Mathematics and Statistics, Huazhong University of Science and Technology,
Wuhan, 430074, China; Hubei Key Laboratory of Engineering Modeling and
Scientific Computing, Huazhong University of Science and Technology, Wuhan
430074, China<EMAIL_ADDRESS>
###### Abstract

In this paper, we first present a multiple-relaxation-time lattice Boltzmann
(MRT-LB) model for one-dimensional diffusion equation where the D1Q3 (three
discrete velocities in one-dimensional space) lattice structure is considered.
Then through the theoretical analysis, we derive an explicit four-level
finite-difference scheme from this MRT-LB model. The results show that the
four-level finite-difference scheme is unconditionally stable, and through
adjusting the weight coefficient $\omega_{0}$ and the relaxation parameters
$s_{1}$ and $s_{2}$ corresponding to the first and second moments, it can also
have a sixth-order accuracy in space. Finally, we also test the four-level
finite-difference scheme through some numerical simulations, and find that the
numerical results are consistent with our theoretical analysis.

###### keywords:

lattice Boltzmann model, diffusion equation, finite-difference scheme.

76M28, 65M06, 65N06, 82C40

##  1 Introduction

In the past decades, the lattice Boltzmann (LB) method, as one of the
mesoscopic numerical approaches, not only has gained a great success in the
study of complex flows governed by the Navier-Stokes equations [1, 2, 3, 4, 5,
6, 7, 8, 9, 10, 11, 12, 13], but also can be considered as a numerical solver
to some other types of the partial differential equations, for example, the
diffusion equations [14, 15, 16, 17, 18], convection-diffusion equations [19,
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], Poisson equation [31, 32, 33],
Kuramoto-Sivashinsky equations [34, 35, 36] and some complex equations [37,
38, 39, 40]. Based on the collision term in the LB method, the LB models can
be classified into three basic kinds, i.e., the lattice BGK or single-
relaxation-time LB (SRT-LB) model [41], two-relaxation-time LB (TRT-LB) model
[20] and the multiple-relaxation-time LB (MRT-LB) model [42]. Actually, the
SRT-LB and TRT-LB models are two special cases of the MRT-LB model [43], and
moreover, the MRT-LB model could be more stable and/or more accurate than SRT-
LB and TRT-LB models through adjusting some free relaxation parameters [44,
45, 46, 47]. For these reasons, the MRT-LB model is considered in this work.

To build the relation between the LB method and the macroscopic partial
differential equations, some asymptotic analysis approaches are usually
adopted [43], including the Chapman-Enskog expansion [48], Maxwell iteration
[49, 50], direct Taylor expansion [51, 52, 53], recurrence equation [56, 57],
and asymptotic expansion with diffusive scaling [58]. With the help of these
asymptotic analysis methods, one can determine the expression of macroscopic
transport efficient, which is related to the relaxation parameter appeared in
the LB model. Although the asymptotic analysis methods can be used to
illustrate that the LB method is suitable for the macroscopic partial
differential equations, and demonstrate the accuracy of the LB method [53, 54,
55], the relation between the mesoscopic LB method and the macroscopic
partial-differential-equation based numerical schemes (hereafter named
macroscopic numerical schemes) is still unclear, and furthermore, can we
obtain the macroscopic numerical scheme of the LB method? Actually, once the
macroscopic numerical scheme of the LB method is obtained, we can not only
gain a better understanding on the LB method through the knowledge already
available on the macroscopic numerical scheme, but also perform more further
research on constructing the mesoscopic LB models for macroscopic partial
differential equations.

In the past years, some efforts have been made on this aspect. Junk [59] and
Inamuro [60] found that when the relaxation parameter is equal to unity, the
SRT-LB model would reduce to a special macroscopic two-level finite-difference
(FD) scheme, and at the diffusive scaling ($\Delta t\propto\Delta x^{2}$,
$\Delta t$ and $\Delta x$ are time step and lattice spacing), the macroscopic
numerical scheme for the incompressible Naiver-Stokes equations has a second-
order convergence rate in space [59]. Ancona [61] demonstrated that for one-
dimensional (1D) convection-diffusion equations, the LB method with D1Q2
lattice structure can be written as the classical DuFort-Frankel scheme [62],
which is an explicit three-level second-order FD scheme. He et al. [63]
performed a theoretical analysis on the SRT-LB model for several simple flows,
and found that when the flows are assumed to be unidirectional and steady-
state, the SRT-LB mode is nothing but a macroscopic second-order FD scheme of
simplified incompressible Navier-Stokes equations. Then they also obtained the
analytical solutions of these simple flows under some commonly used schemes
for nonslip boundary conditions, and demonstrated that usually there is a
numerical slip velocity at the solid wall which cannot be eliminated
effectively in the SRT-LB model. Following the same way, Guo and his
collaborators [64, 65] conducted some analyses on the MRT-LB model for
microscale flows, and obtained the similar results as those reported in Ref.
[63], while the numerical slip caused by the bounce-back scheme can be
eliminated through adjusting the free relaxation parameter corresponding the
third-order moment of the distribution function. d’Humière and Ginzburg [56,
57] also carried out a theoretical analysis on the TRT-LB model with
recurrence equations, and found that when the relation $\Lambda^{eo}=1/4$
($\Lambda^{eo}$ is the so-called _magic_ parameter) is satisfied, the TRT-LB
model would become a macroscopic three-level FD scheme with a second-order
accuracy in space. Li et al. [66] conducted an analysis on the SRT-LB model
with D1Q2 lattice structure for 1D Burgers equation, and demonstrated that
this LB model is equivalent to an explicit three-level second-order FD scheme.
Recently, Cui et al. [46] also carried out a theoretical analysis on the MRT-
LB model for the convection-diffusion equation, and found that for the the
unidirectional and steady-state problems, the MRT-LB model is just a
macroscopic second-order FD scheme of convection-diffusion equation. From the
works mentioned above, one can see that under some special conditions, the LB
method is equivalent to a special macroscopic second-order FD scheme for a
specified partial differential equation, which is also consistent with the
results of some asymptotic analyses on the LB method [43]. However, through
choosing the weight coefficients and relaxation parameter properly, Suga [17]
found that the SRT-LB model with D1Q3 lattice structure could be a macroscopic
four-level fourth-order FD scheme for 1D diffusion equation. We noted that
this work is only limited to the SRT-LB model, it is still unclear whether the
more general MRT-LB model can be written as a macroscopic high-order FD
scheme, and additionally, can we obtain a more accurate FD scheme from the
MRT-LB model where more degrees of freedom in adjusting the relaxation
parameters are included? To answer these questions, in this work we first
develop a MRT-LB model for 1D diffusion equation where a free weight
coefficient $\omega_{0}$ is introduced, and then based on this MRT-LB model,
we obtain an equivalent macroscopic four-level FD scheme. Through some
theoretical analyses, we show that the four-level FD scheme is unconditionally
stable, and can achieve a sixth-order convergence rate in space.

The rest of the paper is organized as follows. In section 2, we presented a
MRT-LB model for the 1D diffusion equations where the D1Q3 lattice structure
is adopted, and then derived an explicit four-level FD scheme from the MRT LB
model. Additionally, it is also shown that the macroscopic numerical schemes
of the SRT-LB model, TRT-LB model, regularized-LB model and modified-lattice-
kinetic model are just some special cases of that of the MRT-LB model. In
section 3, we investigated the accuracy of the four-level FD scheme, followed
by a stability analysis in section 4. In section 5, we performed some
simulations, and found that under some special conditions, the four-level FD
scheme has a sixth-order convergence rate in space, which is also consistent
with our theoretical analysis. Finally, some conclusions are given in section
6.

##  2 The multiple-relaxation-time lattice Boltzmann model based four-level
finite-difference scheme for one-dimensional diffusion equation

In this section, we first developed a MRT-LB model for 1D diffusion equation
with a constant diffusion coefficient, and then presented the details on how
to obtain an explicit four-level FD scheme from the MRT-LB model.

###  2.1 The multiple-relaxation-time lattice Boltzmann model for one-
dimensional diffusion equation

From the mathematical point of view, the 1D diffusion process of mass and heat
can be described by the classical diffusion equation,

(1) $\frac{\partial\phi}{\partial t}=\kappa\frac{\partial^{2}\phi}{\partial
x^{2}}+R,$

where $\phi$ is a scalar variable dependent on the space $x$ and time $t$.
$\kappa$ is the diffusion coefficient, $R$ is the source term, and in this
work, they are assumed to be two constants. In the framework of LB method, the
diffusion equation (1) can be solved efficiently and accurately by different
LB models [14, 15, 17, 20, 25, 28], here we only consider the more general
MRT-LB model for its accuracy and stability in the study of complex problems
[45, 47, 67].

The evolution of MRT-LB model for the diffusion equation (1) can be written as
[28, 46]

(2) $\displaystyle f_{i}(x+\mathbf{c}_{i}\Delta t,t+\Delta t)$
$\displaystyle=f_{i}(x,t)-(\mathbf{M}^{-1}\mathbf{SM})_{ik}[f_{k}(x,t)-f_{k}^{eq}(x,t)]$
$\displaystyle+\Delta
t\big{[}\mathbf{M}^{-1}\big{(}\mathbf{I}-\frac{\mathbf{S}}{2}\big{)}\mathbf{M}\big{]}_{ik}R_{k},(i=-1,0,1),$

where $f_{i}(x,t)$ and $f_{i}^{eq}(x,t)$ are the distribution function and
equilibrium distribution at position $x$ and time $t$. In the D1Q3 lattice
structure, the discrete velocity $\mathbf{c}_{i}$, the transformation matrix
$\mathbf{M}$ and the diagonal relaxation matrix $\mathbf{S}$ can be given by

(3a) $\mathbf{c}_{i}=\left\\{\begin{array}[]{ll}-c,&\mathrm{\emph{i}=-1},\\\
0,&\mathrm{\emph{i}=0},\\\ c,&\mathrm{\emph{i}=1},\end{array}\right.$ (3b)
$\mathbf{M}=\left(\begin{array}[]{ccc}1&1&1\\\ -c&0&c\\\
c^{2}&-2c^{2}&c^{2}\end{array}\right),$ (3c)
$\mathbf{S}=\left(\begin{array}[]{ccc}s_{0}&0&0\\\ 0&s_{1}&0\\\
0&0&s_{2}\end{array}\right),$

where $c=\Delta x/\Delta t$ is the lattice speed with $\Delta x$ and $\Delta
t$ being lattice spacing and time step, respectively. The diagonal element
$s_{i}$ of the relaxation matrix $\mathbf{S}$ is the relaxation parameter
corresponding to $i$th moment of the distribution function $f_{i}(x,t)$, and
to make the physical transport coefficient (e.g., diffusion coefficient)
positive, it should be located in the range $(0,\ 2)$. $R_{i}$ is the discrete
source term, and is defined by

(4) $R_{i}=\omega_{i}R,$

where $\omega_{i}$ is the weight coefficient.

In the LB method, to derive correct macroscopic diffusion equation (1), the
equilibrium distribution should be defined as

(5) $f_{i}^{eq}(x,t)=\omega_{i}\phi(x,t),$

which satisfies the following conditions [28],

(6a) $\sum_{i}f_{i}^{eq}(x,\ t)=\phi(x,\ t),$ (6b)
$\sum_{i}\mathbf{c}_{i}f_{i}^{eq}(x,\ t)=0,$ (6c)
$\sum_{i}\mathbf{c}_{i}\mathbf{c}_{i}f_{i}^{eq}(x,\ t)=2\omega_{1}\phi(x,\
t)c^{2},$

where the relation $\omega_{-1}=\omega_{1}$ derived from Eq. (6b) is used to
obtain Eq. (6c). If the weight coefficient $\omega_{0}$ is considered as a
free parameter, then from Eq. (6a) and the condition $\omega_{-1}=\omega_{1}$,
one can obtain

(7) $\omega_{1}=\omega_{-1}=\frac{1-\omega_{0}}{2},$

where $0<\omega_{0}<1$ which can be used to ensure that all weight
coefficients are larger than zero. In addition, the macroscopic variable
$\phi(x,\ t)$ can be calculated by

(8) $\phi(x,t)=\sum_{i}f_{i}(x,\ t)+\frac{\Delta t}{2}R.$

Through the Chapman-Enskog analysis [28], one can correctly recover the
diffusion equation (1) from the present MRT-LB model with the following
relation between the diffusion coefficient and relaxation parameter $s_{1}$,

(9) $\kappa=2\omega_{1}\big{(}\frac{1}{s_{1}}-\frac{1}{2}\big{)}\frac{\Delta
x^{2}}{\Delta t}.$

###  2.2 The multiple-relaxation-time lattice Boltzmann model based explicit
four-level finite-difference scheme

In this part, we will show some details on how to derive the macroscopic
numerical scheme from the present MRT-LB model. To simplify the following
analysis, the notations $f_{i,j}^{n}=f_{i}(j\Delta x,\ n\Delta t)$ and
$\phi_{j}^{n}=\phi(j\Delta x,\ n\Delta t)$ are introduced. Through
substituting the discrete source term $R_{i}$ [Eq. (4)] and equilibrium
distribution function [Eq. (5)] into the evolution equation (2), we have

(10a)
$f_{-1,j}^{n+1}=f_{-1,j+1}^{n}-\frac{s_{1}}{2}\big{(}f_{-1,j+1}^{n}-f_{1,j+1}^{n}\big{)}+\frac{s_{2}}{2}f_{0,j+1}^{n}-\frac{\omega_{0}s_{2}}{2}\phi_{j+1}^{n}+\big{(}\omega_{-1}+\frac{\omega_{0}s_{2}}{4}\big{)}\Delta
tR,$ (10b)
$f_{0,j}^{n+1}=f_{0,j}^{n}-s_{2}f_{0,j}^{n}+\omega_{0}s_{2}\phi_{j}^{n}+\omega_{0}\big{(}1-\frac{s_{2}}{2}\big{)}\Delta
tR,$ (10c)
$f_{1,j}^{n+1}=f_{1,j-1}^{n}+\frac{s_{1}}{2}\big{(}f_{-1,j-1}^{n}-f_{1,j-1}^{n}\big{)}+\frac{s_{2}}{2}f_{0,j-1}^{n}-\frac{\omega_{0}s_{2}}{2}\phi_{j-1}^{n}+\big{(}\omega_{1}+\frac{\omega_{0}s_{2}}{4}\big{)}\Delta
tR,$

where Eq. (8) has been used. To obtain the macroscopic numerical scheme of the
MRT-LB model, the distribution functions appeared in Eq. (10) must be replaced
by the macroscopic variable $\phi$ at different grid points and time levels.
For this purpose, we first conduct a sum of Eq. (10), and derive the following
equation,

(11) $\displaystyle\phi_{j}^{n+1}$
$\displaystyle=f_{-1,j+1}^{n}+f_{{0,j}}^{n}+f_{1,j-1}^{n}-\frac{s_{1}}{2}\big{[}\big{(}f_{-1,j+1}^{n}-f_{1,j+1}^{n}\big{)}-\big{(}f_{-1,j-1}^{n}-f_{1,j-1}^{n}\big{)}\big{]}$
$\displaystyle+\frac{s_{2}}{2}\big{(}f_{0,j-1}^{n}-2f_{0,j}^{n}+f_{0,j+1}^{n}\big{)}-\frac{\omega_{0}s_{2}}{2}\big{(}\phi_{j+1}^{n}-2\phi_{j}^{n}+\phi_{j-1}^{n}\big{)}+\frac{3}{2}\Delta
tR.$

Then from Eq. (10b) we can obtain

(12)
$-f_{0,j}^{n}+s_{2}f_{0,j}^{n}-\omega_{0}s_{2}\phi_{j}^{n}=-f_{0,j}^{n+1}+\omega_{0}\big{(}1-\frac{s_{2}}{4}\big{)}\Delta
tR.$

Based on Eqs. (8), (10a), (10c) and (12), we have

(13)
$f_{1,j+1}^{n}+f_{-1,j-1}^{n}=-f_{0,j}^{n}+\phi_{j}^{n-1}+\frac{1}{2}\Delta
tR.$

In addition, from (8) one can derive

(14) $\displaystyle f_{-1,j+1}^{n}+f_{0,j}^{n}+f_{1,j-1}^{n}$
$\displaystyle=-f_{-1,j-1}^{n}-f_{1,j+1}^{n}-f_{0,j-1}^{n}-f_{0,j+1}^{n}+f_{0,j}^{n}$
$\displaystyle+\phi_{j+1}^{n}+\phi_{j-1}^{n}-\Delta tR$
$\displaystyle=-\big{(}f_{0,j-1}^{n}-2f_{0,j}^{n}+f_{0,j+1}^{n}\big{)}+\phi_{j+1}^{n}+\phi_{j-1}^{n}-\phi_{j}^{n-1}-\frac{3}{2}\Delta
tR,$

where Eq. (13) has been adopted. Similarly, from Eqs. (8) and (13) we can also
obtain

(15)
$\displaystyle\big{(}f_{-1,j+1}^{n}-f_{1,j+1}^{n}\big{)}-\big{(}f_{-1,j-1}^{n}-f_{1,j-1}^{n}\big{)}$
$\displaystyle=-2\big{(}f_{1,j+1}^{n}+f_{-1,j-1}^{n}\big{)}-f_{0,j+1}^{n}-f_{0,j-1}^{n}+\phi_{j+1}^{n}+\phi_{j-1}^{n}-\Delta
tR$
$\displaystyle=-\big{(}f_{0,j-1}^{n}-2f_{0,j}^{n}+f_{0,j+1}^{n}\big{)}+\phi_{j+1}^{n}+\phi_{j-1}^{n}-2\phi_{j}^{n-1}-2\Delta
tR.$

Substituting Eqs. (14) and (LABEL:eq2-15) into Eq. (11) yields

(16) $\displaystyle\phi_{j}^{n+1}=$
$\displaystyle\big{(}\frac{s_{1}}{2}+\frac{s_{2}}{2}-1\big{)}\big{(}f_{0,j-1}^{n}-2f_{0,j}^{n}+f_{0,j+1}^{n}\big{)}+\big{(}1-\frac{s_{1}}{2}-\frac{\omega_{0}s_{2}}{2}\big{)}\phi_{j+1}^{n}$
$\displaystyle+\omega_{0}s_{2}\phi_{j}^{n}+\big{(}1-\frac{s_{1}}{2}-\frac{\omega_{0}s_{2}}{2}\big{)}\phi_{j-1}^{n}+(s_{1}-1)\phi_{j}^{n-1}+\Delta
ts_{1}R.$

Now we need to give an evaluation of the first term on the right hand side of
Eq. (16). Actually, from Eqs. (8) and (10) we have

(17) $\displaystyle f_{0,j-1}^{n}-2f_{0,j}^{n}+f_{0,j+1}^{n}=$ $\displaystyle
f_{0,j-1}^{n}+f_{0,j+1}^{n}-2(\phi_{j}^{n}-f_{-1,j}^{n}-f_{1,j}^{n})+\Delta
tR$ $\displaystyle=$ $\displaystyle
2(f_{-1,j}^{n}+f_{1,j}^{n})+f_{0,j-1}^{n}+f_{0,j+1}^{n}-2\phi_{j}^{n}+\Delta
tR$ $\displaystyle=$ $\displaystyle
f_{0,j+1}^{n-1}+2f_{-1,j+1}^{n-1}+f_{0,j-1}^{n-1}+2f_{1,j-1}^{n-1}$
$\displaystyle-$ $\displaystyle
s_{1}\big{[}\big{(}f_{-1,j+1}^{n-1}-f_{1,j+1}^{n-1}\big{)}-\big{(}f_{-1,j-1}^{n-1}-f_{1,j-1}^{n-1}\big{)}\big{]}-2\phi_{j}^{n}+3\Delta
tR.$

With the help of Eqs. (8) and (13), one can derive the following equation,

(18) $\displaystyle f_{0,j+1}^{n-1}$
$\displaystyle+2f_{-1,j+1}^{n-1}+f_{0,j-1}^{n-1}+2f_{1,j-1}^{n-1}$
$\displaystyle=$ $\displaystyle-
f_{0,j+1}^{n-1}-f_{0,j-1}^{n-1}-2\big{(}f_{-1,j}^{n-1}+f_{1,j}^{n-1}\big{)}+2\phi_{j+1}^{n-1}+2\phi_{j-1}^{n-1}-2\Delta
tR$ $\displaystyle=$
$\displaystyle-\big{(}f_{0,j-1}^{n-1}-2f_{0,j}^{n-1}+f_{0,j+1}^{n-1}\big{)}+2\phi_{j+1}^{n-1}+2\phi_{j-1}^{n-1}-2\phi_{j}^{n-2}-3\Delta
tR.$

Substituting Eqs. (LABEL:eq2-15) and (18) into Eq. (17) gives rise to

(19) $\displaystyle f_{0,j-1}^{n}-2f_{0,j}^{n}+f_{0,j+1}^{n}=$
$\displaystyle(s_{1}-1)\big{(}f_{0,j-1}^{n-1}-2f_{0,j}^{n-1}+f_{0,j+1}^{n-1}\big{)}$
$\displaystyle-$ $\displaystyle
2\phi_{j}^{n}+(2-s_{1})\phi_{j-1}^{n-1}+(2-s_{1})\phi_{j}^{n-2}-(2s_{1}-2)\phi_{j}^{n-2}$
$\displaystyle+$ $\displaystyle 2s_{1}\Delta tR.$

If we substitute Eq. (19) into Eq. (16), one can obtain

(20) $\displaystyle\phi_{j}^{n+1}=$
$\displaystyle(s_{1}-1)\big{(}\frac{s_{1}}{2}+\frac{s_{2}}{2}-1\big{)}\big{(}f_{{0}_{j-1}}^{n-1}-2f_{{0}_{j}}^{n-1}+f_{{0}_{j+1}}^{n-1}\big{)}$
$\displaystyle+$
$\displaystyle\big{(}1-\frac{s_{1}}{2}-\frac{\omega_{0}s_{2}}{2}\big{)}\phi_{j-1}^{n}+[\omega_{0}s_{2}-2\big{(}\frac{s_{1}}{2}+\frac{s_{2}}{2}-1\big{)}\big{]}\phi_{j}^{n}$
$\displaystyle+$
$\displaystyle\big{(}1-\frac{s_{1}}{2}-\frac{\omega_{0}s_{2}}{2}\big{)}\phi_{j+1}^{n}+(2-s_{1})\big{(}\frac{s_{1}}{2}+\frac{s_{2}}{2}-1\big{)}\phi_{j-1}^{n-1}+(s_{1}-1)\phi_{j}^{n-1}$
$\displaystyle+$
$\displaystyle(2-s_{1})\big{(}\frac{s_{1}}{2}+\frac{s_{2}}{2}-1\big{)}\phi_{j+1}^{n-1}+(2s_{1}-2)\big{(}\frac{s_{1}}{2}+\frac{s_{2}}{2}-1\big{)}\phi_{j}^{n-2}$
$\displaystyle+$ $\displaystyle s_{1}(s_{1}+s_{2}-1)\Delta tR.$

With the aid of Eq. (16), we can rewrite Eq. (20) as

(21) $\displaystyle\phi_{j}^{n+1}=$
$\displaystyle\big{(}1-\frac{s_{1}}{2}-\frac{\omega_{0}s_{2}}{2}\big{)}\phi_{j-1}^{n}+[(\omega_{0}-1)s_{2}+1]\phi_{j}^{n}+\big{(}1-\frac{s_{1}}{2}-\frac{\omega_{0}s_{2}}{2}\big{)}\phi_{j+1}^{n}$
$\displaystyle+$
$\displaystyle\big{(}\frac{\omega_{0}s_{1}s_{2}}{2}-\frac{s_{1}s_{2}}{2}-\frac{\omega_{0}s_{2}}{2}+\frac{s_{1}}{2}+s_{2}-1\big{)}\phi_{j-1}^{n-1}$
$\displaystyle+$
$\displaystyle(-\omega_{0}s_{1}s_{2}+\omega_{0}s_{2}+s_{1}-1)\phi_{j}^{n-1}$
$\displaystyle+$
$\displaystyle\big{(}\frac{\omega_{0}s_{1}s_{2}}{2}-\frac{s_{1}s_{2}}{2}-\frac{\omega_{0}s_{2}}{2}+\frac{s_{1}}{2}+s_{2}-1\big{)}\phi_{j+1}^{n-1}$
$\displaystyle+$
$\displaystyle(s_{1}-1)(s_{2}-1)\phi_{j}^{n-2}+s_{1}s_{2}\Delta tR,$

which can also be written as

(22) $\displaystyle\phi_{j}^{n+1}=$
$\displaystyle\alpha_{1}\phi_{j-1}^{n}+\alpha_{2}\phi_{j}^{n}+\alpha_{1}\phi_{j+1}^{n}$
$\displaystyle+$
$\displaystyle\beta_{1}\phi_{j-1}^{n-1}+\beta_{2}\phi_{j}^{n-1}+\beta_{1}\phi_{j+1}^{n-1}$
$\displaystyle+$ $\displaystyle\gamma\phi_{j}^{n-2}+\delta\Delta tR,$

where the parameters $\alpha_{i}$ ($i=1,2$), $\beta_{i}$, $\gamma$ and
$\delta$ are given by

(23) $\displaystyle\alpha_{1}$
$\displaystyle=1-\frac{s_{1}}{2}-\frac{\omega_{0}s_{2}}{2},$
$\displaystyle\alpha_{2}$ $\displaystyle=(\omega_{0}-1)s_{2}+1,$
$\displaystyle\beta_{1}$
$\displaystyle=\frac{\omega_{0}s_{1}s_{2}}{2}-\frac{s_{1}s_{2}}{2}-\frac{\omega_{0}s_{2}}{2}+\frac{s_{1}}{2}+s_{2}-1,$
$\displaystyle\beta_{2}$
$\displaystyle=-\omega_{0}s_{1}s_{2}+\omega_{0}s_{2}+s_{1}-1,$
$\displaystyle\gamma$ $\displaystyle=(s_{1}-1)(s_{2}-1),$
$\displaystyle\delta$ $\displaystyle=s_{1}s_{2}.$

Here we would like to point out that Eq. (22) is the the mesoscopic MRT-LB
model based macroscopic explicit four-level FD scheme for 1D diffusion
equation.
Remark I: If $s_{0}=s_{1}=s_{2}=\omega$ ($\omega$ is the relaxation parameter
in the SRT-LB model) and $R=0$, we can obtain the following SRT-LB model based
four-level FD scheme from Eq. (22),

(24) $\displaystyle\phi_{j}^{n+1}=$
$\displaystyle\alpha_{1}\phi_{j-1}^{n}+\alpha_{2}\phi_{j}^{n}+\alpha_{1}\phi_{j+1}^{n}$
$\displaystyle+$
$\displaystyle\beta_{1}\phi_{j-1}^{n-1}+\beta_{2}\phi_{j}^{n-1}+\beta_{1}\phi_{j+1}^{n-1}$
$\displaystyle+$ $\displaystyle\gamma\phi_{j}^{n-2},$

with the following parameters,

(25) $\displaystyle\alpha_{1}=\Omega+\omega_{1}\omega,$
$\displaystyle\alpha_{2}=\Omega+(1-2\omega_{1})\omega,$
$\displaystyle\beta_{1}=-[\Omega+(1-\omega_{1})\omega]\Omega,\ \ \beta_{2}$
$\displaystyle=-(\Omega+2\omega_{1}\omega)\Omega,\ \ \gamma$
$\displaystyle=\Omega^{2},$

where $\Omega=1-\omega$. It is clear that the SRT-LB model based macroscopic
numerical scheme (24), as a special case of the present MRT-LB model based
four-level scheme, is the same as that reported in the previous work [17].
Remark II: If $s_{0}=s_{2}=s^{+}$ and $s_{1}=s^{-}$ with $s^{+}$ and $s^{-}$
denoting the relaxation parameters corresponding to the symmetric and
antisymmetric modes [20], we can obtain the TRT-LB model based four-level FD
scheme from Eq. (22). However, when $s_{0}=s_{2}=1$ and $s_{1}=\omega$, one
can derive the regularized-LB model [68, 69] based macroscopic numerical
scheme from Eq. (22). In addition, if $s_{0}=s_{2}=\omega$ and
$s_{1}=\omega/(1-\omega\eta)$ with $\eta$ being an adjusting parameter, we can
also obtain the modified-lattice-kinetic model [70, 71, 72] based macroscopic
numerical scheme from Eq. (22).

##  3 The accuracy analysis of the multiple-relaxation-time lattice Boltzmann
model based macroscopic numerical scheme

We now performed an accuracy analysis on the macroscopic four-level FD scheme
(22). To do this, we first conducted the Taylor expansion to Eq. (22) at the
position $x=j\Delta x$ and time $t=n\Delta t$, and after some algebraic
manipulations, one can obtained

(26)
$\displaystyle(1+2\beta_{1}+\beta_{2}+2\gamma)\big{[}\frac{\partial\phi}{\partial
t}\big{]}_{j}^{n}=$ $\displaystyle(\alpha_{1}+\beta_{1})\frac{\Delta
x^{2}}{\Delta t}\big{[}\frac{\partial^{2}\phi}{\partial
x^{2}}\big{]}_{j}^{n}+\delta R$ $\displaystyle+$
$\displaystyle\frac{1}{12}(\alpha_{1}+\beta_{1})\frac{\Delta x^{4}}{\Delta
t}\big{[}\frac{\partial^{4}\phi}{\partial
x^{4}}\big{]}_{j}^{n}+\frac{1}{360}(\alpha_{1}+\beta_{1})\frac{\Delta
x^{6}}{\Delta t}\big{[}\frac{\partial^{6}\phi}{\partial x^{6}}\big{]}_{j}^{n}$
$\displaystyle-$ $\displaystyle\beta_{1}\Delta
x^{2}\big{[}\frac{\partial^{3}\phi}{\partial x^{2}\partial
t}\big{]}_{j}^{n}-\frac{\beta_{1}}{12}\Delta
x^{4}\big{[}\frac{\partial^{5}\phi}{\partial x^{4}\partial t}\big{]}_{j}^{n}$
$\displaystyle+$ $\displaystyle\frac{\beta_{1}}{2}\Delta x^{2}\Delta
t\big{[}\frac{\partial^{4}\phi}{\partial x^{2}\partial
t^{2}}\big{]}_{j}^{n}+\frac{1}{2}\big{(}2\beta_{1}+\beta_{2}+4\gamma-1\big{)}\Delta
t\big{[}\frac{\partial^{2}\phi}{\partial t^{2}}\big{]}_{j}^{n}$
$\displaystyle-$
$\displaystyle\frac{1}{6}\big{(}2\beta_{1}+\beta_{2}+8\gamma+1\big{)}\Delta
t^{2}\big{[}\frac{\partial^{3}\phi}{\partial t^{3}}\big{]}_{j}^{n}+\cdots.$

Substituting Eq. (23) into above equation yields

(27) $\displaystyle\big{[}\frac{\partial\phi}{\partial t}\big{]}_{j}^{n}=$
$\displaystyle\kappa\big{[}\frac{\partial^{2}\phi}{\partial
x^{2}}\big{]}_{j}^{n}+R$ $\displaystyle+$
$\displaystyle\frac{1}{12}\kappa\Delta
x^{2}\big{[}\frac{\partial^{4}\phi}{\partial
x^{4}}\big{]}_{j}^{n}+\frac{1}{360}\kappa\Delta
x^{4}\big{[}\frac{\partial^{6}\phi}{\partial x^{6}}\big{]}_{j}^{n}$
$\displaystyle-$
$\displaystyle\big{(}\frac{\omega_{0}}{2}-\frac{1}{2}-\frac{\omega_{0}}{2s_{1}}+\frac{1}{2s_{2}}+\frac{1}{s_{1}}-\frac{1}{s_{1}s_{2}}\big{)}\Delta
x^{2}\big{[}\frac{\partial^{3}\phi}{\partial x^{2}\partial t}\big{]}_{j}^{n}$
$\displaystyle-$
$\displaystyle\big{(}\frac{\omega_{0}}{24}-\frac{1}{24}-\frac{\omega_{0}}{24s_{1}}+\frac{1}{24s_{2}}+\frac{1}{12s_{1}}-\frac{1}{12s_{1}s_{2}}\big{)}\Delta
x^{4}\big{[}\frac{\partial^{5}\phi}{\partial x^{4}\partial t}\big{]}_{j}^{n}$
$\displaystyle+$
$\displaystyle\big{(}\frac{\omega_{0}}{4}-\frac{1}{4}-\frac{\omega_{0}}{4s_{1}}+\frac{1}{4s_{2}}+\frac{1}{2s_{1}}-\frac{1}{2s_{1}s_{2}}\big{)}\Delta
x^{2}\Delta t\big{[}\frac{\partial^{4}\phi}{\partial x^{2}\partial
t^{2}}\big{]}_{j}^{n}$ $\displaystyle+$
$\displaystyle\big{(}\frac{3}{2}-\frac{1}{s_{1}}-\frac{1}{s_{2}}\big{)}\Delta
t\big{[}\frac{\partial^{2}\phi}{\partial
t^{2}}\big{]}_{j}^{n}+\big{(}\frac{1}{s_{2}}+\frac{1}{s_{1}}-\frac{7}{6}-\frac{1}{s_{1}s_{2}}\big{)}\Delta
t^{2}\big{[}\frac{\partial^{3}\phi}{\partial t^{3}}\big{]}_{j}^{n}+\cdots,$

where Eq. (9) has been used to derive above equation.

According to the diffusion equation (1), we have the following relations,

(28) $\displaystyle\big{[}\frac{\partial^{2}\phi}{\partial
t^{2}}\big{]}_{j}^{n}=\kappa^{2}\big{[}\frac{\partial^{4}\phi}{\partial
x^{4}}\big{]}_{j}^{n},$ $\displaystyle\big{[}\frac{\partial^{3}\phi}{\partial
t^{3}}\big{]}_{j}^{n}=\kappa^{3}\big{[}\frac{\partial^{6}\phi}{\partial
x^{6}}\big{]}_{j}^{n},$ $\displaystyle\big{[}\frac{\partial^{5}\phi}{\partial
x^{4}\partial t}\big{]}_{j}^{n}=\kappa\big{[}\frac{\partial^{6}\phi}{\partial
x^{6}}\big{]}_{j}^{n},\ \ \big{[}\frac{\partial^{3}\phi}{\partial
x^{2}\partial t}\big{]}_{j}^{n}$
$\displaystyle=\kappa\big{[}\frac{\partial^{4}\phi}{\partial
x^{4}}\big{]}_{j}^{n},\ \ \big{[}\frac{\partial^{4}\phi}{\partial
x^{2}\partial
t^{2}}\big{]}_{j}^{n}=\kappa^{2}\big{[}\frac{\partial^{6}\phi}{\partial
x^{6}}\big{]}_{j}^{n},$

Substituting Eq. (28) into Eq. (27), we have

(29) $\displaystyle\big{[}\frac{\partial\phi}{\partial t}\big{]}_{j}^{n}=$
$\displaystyle\kappa\big{[}\frac{\partial^{2}\phi}{\partial
x^{2}}\big{]}_{j}^{n}+R$ $\displaystyle+$
$\displaystyle\frac{1}{s_{1}s_{2}}\big{[}\frac{s_{1}s_{2}}{12}-\big{(}\frac{\omega_{0}s_{1}s_{2}}{2}-\frac{s_{1}s_{2}}{2}-\frac{\omega_{0}s_{2}}{2}+\frac{s_{1}}{2}+s_{2}-1\big{)}$
$\displaystyle+$
$\displaystyle\big{(}\frac{3s_{1}s_{2}}{2}-s_{2}-s_{1}\big{)}\epsilon\big{]}\kappa\Delta
x^{2}\big{[}\frac{\partial^{4}\phi}{\partial x^{4}}\big{]}_{j}^{n}$
$\displaystyle+$
$\displaystyle\frac{1}{s_{1}s_{2}}\big{[}\frac{s_{1}s_{2}}{360}-\big{(}\frac{\omega_{0}s_{1}s_{2}}{24}-\frac{s_{1}s_{2}}{24}-\frac{\omega_{0}s_{2}}{24}+\frac{s_{1}}{24}+\frac{s_{2}}{12}-\frac{1}{12}\big{)}$
$\displaystyle+$
$\displaystyle\big{(}\frac{\omega_{0}s_{1}s_{2}}{4}-\frac{s_{1}s_{2}}{4}-\frac{\omega_{0}s_{2}}{4}+\frac{s_{1}}{4}+\frac{s_{2}}{2}-\frac{1}{2}\big{)}\epsilon$
$\displaystyle+$
$\displaystyle\big{(}-\frac{7s_{1}s_{2}}{6}+s_{2}+s_{1}-1\big{)}\epsilon^{2}\big{]}\kappa\Delta
x^{4}\big{[}\frac{\partial^{6}\phi}{\partial x^{6}}\big{]}_{j}^{n}$
$\displaystyle+$ $\displaystyle O(\Delta x^{6}+\Delta t^{3}),$

where $\epsilon$ is the mesh Fourier number, and is defined by

(30) $\epsilon=\frac{\kappa\Delta t}{\Delta
x^{2}}=2\omega_{1}\big{(}\frac{1}{s_{1}}-\frac{1}{2}\big{)}=(1-\omega_{0})\big{(}\frac{1}{s_{1}}-\frac{1}{2}\big{)}.$

From above discussion, it is clear that for a given diffusion coefficient
$\kappa$ or mesh Fourier number $\epsilon$, one can obtain an explicit four-
level FD scheme with the third-order accuracy in time and sixth-order accuracy
in space once the following conditions are satisfied, i.e., the second and
fourth-order truncation errors in Eq. (29) are equal to zero,

(31a)
$\frac{s_{1}s_{2}}{12}-\big{(}\frac{\omega_{0}s_{2}}{2}+\frac{s_{1}}{2}-1\big{)}+\big{(}\frac{s_{1}s_{2}}{2}-s_{2}-s_{1}\big{)}\epsilon=0\\\
$ (31b)
$\frac{s_{1}s_{2}}{360}-\frac{1}{12}\big{(}\frac{\omega_{0}s_{2}}{2}+\frac{s_{1}}{2}-1\big{)}-\frac{1}{2}\big{(}\frac{s_{1}s_{2}}{6}-\frac{\omega_{0}s_{2}}{2}-\frac{s_{1}}{2}+1\big{)}\epsilon+\big{(}-\frac{2s_{1}s_{2}}{3}+s_{2}+s_{1}-1\big{)}\epsilon^{2}=0,$

where Eq. (30) has been applied. We note that compared to the SRT-LB model
based four-level FD scheme for 1D equation [17], the present MRT-LB model
based macroscopic numerical scheme can be more accurate through adjusting the
weight coefficient $\omega_{0}$ and relaxation parameters $s_{1}$ and $s_{2}$
to satisfy Eq. (31). In the following, some remarks on the MRT-LB model based
macroscopic numerical scheme are listed.
Remark I: It is found that the relaxation parameter $s_{0}$ corresponding to
the zeroth moment of distribution function $f_{i}(x,t)$ (or conservative
variable $\phi$) does not appear in the macroscopic numerical scheme (22), and
thus it has no influence on the numerical results. This also explains why the
relaxation parameter $s_{0}$ in the MRT-LB method can be chosen arbitrarily.
However, unlike the relaxation parameter $s_{0}$, the relaxation parameter
$s_{2}$ corresponding to the second-order moment of distribution function has
an important influence on the macroscopic numerical scheme (22) [see Eq.
(31)], and also the numerical results. Actually, these results have also been
reported in some previous works on the MRT-LB model for (convection) diffusion
equations [28, 46, 67].
Remark II: When the conditions of Eqs. (30) and (31) are satisfied, the MRT-LB
model based macroscopic four-level numerical scheme (22) has a sixth-order
accuracy in space. However, if only Eqs. (30) and (31a) hold, the macroscopic
four-level numerical scheme would have a fourth-order accuracy in space.
Compared to the SRT-LB model based macroscopic four-level fourth-order
numerical scheme [17], the present macroscopic numerical scheme has another
distinct characteristic, i.e., there is a free relaxation parameter $s_{2}$,
which can be used to remove the discrete effect of anti-bounce-back scheme for
the Dirichlet boundary conditions [46, 67].
Remark III: Compared to the classical two-level LB method, the implementation
of the macroscopic four-level numerical scheme needs the initial values of
variable $\phi$ at first three time levels, and usually to complete the
initialization, we must adopt some other numerical schemes to obtain the
values of variable $\phi$ at the second and third time levels. In addition,
the inclusion of the extra time levels also brings a larger memory requirement
to store the variable $\phi$.

Finally, to implement the MRT-LB model based macroscopic sixth-order numerical
scheme (22), we must determine the weight coefficient $\omega_{0}$ and
relaxation parameters $s_{1}$ and $s_{2}$ from Eqs. (30) and (31) for a
specified $\epsilon$, while due to the coupling and nonlinearity of these
equations, it is difficult to derive the explicit expressions of $\omega_{0}$,
$s_{1}$ and $s_{2}$ in terms of $\epsilon$. For this reason, we numerically
solved Eqs. (30) and (31), and plotted $\omega_{0}$, $s_{1}$ and $s_{2}$ as a
function of the mesh Fourier number $\epsilon$ in Fig. 1. Here it should be
noted that only under the condition of $0<\epsilon\leq\epsilon_{max}$ with
$\epsilon_{max}$ being about 0.262, we can obtain the real roots of Eqs. (30)
and (31). As seen from Fig. 1, for a given $\epsilon$, the weight coefficient
$\omega_{0}$ and relaxation parameters $s_{1}$ and $s_{2}$ are located in the
ranges of $0.8\leq\omega_{0}<1$, $0<s_{1}\leq s_{1,max}$ and $s_{2,min}\leq
s_{2}<2.0$ with $s_{1,max}$ and $s_{2,min}$ being very close to 0.920 and
1.124.

Figure 1: The weight coefficient $\omega_{0}$ and the relaxation parameters
$s_{1}$ and $s_{2}$ as a function of the discretization parameter $\epsilon$.

##  4 The stability analysis of the multiple-relaxation-time lattice Boltzmann
model based macroscopic numerical scheme

In this section, we will prove that the MRT-LB model based four-level FD
scheme (22) is unconditionally stable. To this end, we first neglected the
source term $R$ and replaced the $\phi_{j}^{n}$ in Eq. (10) by the
distribution function $f_{i,j}^{n}$ through the relation (8), and then take
the discrete Fourier transform of $f_{i,j}^{n}$ in Eq. (10) to obtain the
following matrix equation [73],

(32)
$\widehat{\mathbf{U}}_{j}^{n+1}=\mathbf{G}(\theta,\omega_{0},s_{1},s_{2})\widehat{\mathbf{U}}_{j}^{n},$

where $\widehat{\mathbf{U}}_{j}^{n}$ is the discrete Fourier transform of
$f_{i,j}^{n}$ $(i=-1,0,1)$. $\mathbf{G}$ is the amplification matrix of the
scheme, and is given by

(33)
$\mathbf{G}=\left(\begin{array}[]{cccc}\big{(}1-\frac{s_{1}}{2}-\frac{\omega_{0}s_{2}}{2}\big{)}e^{i\theta}&\big{(}\frac{s_{2}}{2}-\frac{\omega_{0}s_{2}}{2}\big{)}e^{i\theta}&\big{(}\frac{s_{1}}{2}-\frac{\omega_{0}s_{2}}{2}\big{)}e^{i\theta}\\\
\omega_{0}s_{2}&\omega_{0}s_{2}-s_{2}+1&\omega_{0}s_{2}\\\
\big{(}\frac{s_{1}}{2}-\frac{\omega_{0}s_{2}}{2}\big{)}e^{-i\theta}&\big{(}\frac{s_{2}}{2}-\frac{\omega_{0}s_{2}}{2}\big{)}e^{-i\theta}&\big{(}1-\frac{s_{1}}{2}-\frac{\omega_{0}s_{2}}{2}\big{)}e^{-i\theta}\end{array}\right),$

where $-\pi\leq\theta\leq\pi$. On the other hand, if we take the discrete
Fourier transform of Eq. (22), one can obtain the amplification matrix
$\mathbf{H}$ of the macroscopic numerical scheme,

(34)
$\mathbf{H}=\left(\begin{array}[]{cccc}2\alpha_{1}\cos\theta+\alpha_{2}&2\beta_{1}\cos\theta+\beta_{2}&\gamma\\\
1&0&0\\\ 0&1&0\end{array}\right).$

Although the amplification matrix $\mathbf{H}$ is different from $\mathbf{G}$,
due to the equivalence between the MRT-LB model [Eq. (2) or (10)] and the
macroscopic numerical scheme (22), the characteristic polynomials of them are
identical, and can be expressed as

(35) $p(\lambda)=\lambda^{3}+p_{2}\lambda^{2}+p_{1}\lambda+p_{0},$

where the coefficients $p_{0}$, $p_{1}$ and $p_{2}$ are given by

(36) $\displaystyle p_{0}$ $\displaystyle=(s_{1}-1)(1-s_{2}),$ $\displaystyle
p_{1}$
$\displaystyle=(s_{1}-1)(s_{2}\omega_{0}-1)+[(s_{1}-2)(s_{2}-1)+s_{2}\omega_{0}(1-s_{1})]\cos\theta,$
$\displaystyle p_{2}$
$\displaystyle=s_{2}-s_{2}\omega_{0}-1+(s_{2}\omega_{0}+s_{1}-2)\cos\theta.$

In the following, we would show that the roots of the characteristic
polynomial $p(\lambda)$ denoted by $\lambda_{k}$ ($k$=1, 2 and 3) satisfy the
condition $|\lambda_{k}|\leq 1$.

With the following linear fractional transformation,

(37) $\lambda=\frac{1+z}{1-z},\ z\in\mathcal{C},$

the unit circle $|\lambda|=1$ and the field $|\lambda|<1$ are mapped to the
imaginary axis [Re$(z)=0$] and left-half plane [Re$(z)<0$], and vise verse.
Here $\mathcal{C}$ and Re denote the complex-number field and the real part of
a complex number. Substituting Eq. (37) into Eq. (35), we have

(38) $\displaystyle(1-z)^{3}p\big{(}\frac{1+z}{1-z}\big{)}$
$\displaystyle=(1+z)^{3}+p_{2}(1-z)(1+z)^{2}+p_{1}(1-z)^{2}(1+z)+p_{0}(1-z)^{3}$
$\displaystyle=(1-p_{0}+p_{1}-p_{2})z^{3}+(3+3p_{0}-p_{1}-p_{2})z^{2}$
$\displaystyle+(3-3p_{0}-p_{1}+p_{2})z+(1+p_{0}+p_{1}+p_{2}).$

To ensure that the roots of characteristic polynomial $p(\lambda)$ are located
in the field $|\lambda|<1$, the following conditions must be satisfied, i.e.,
the Routh-Hurwitz stability criterion [74, 75, 76],

(39a) $1-p_{0}+p_{1}-p_{2}>0$ (39b) $3+3p_{0}-p_{1}-p_{2}>0,$ (39c)
$3-3p_{0}-p_{1}+p_{2}>0,$ (39d) $1+p_{0}+p_{1}+p_{2}>0,$ (39e)
$1-p_{1}+p_{0}p_{2}-p_{0}^{2}>0.$

After the sums of Eqs. (39a) and (39c), and Eqs. (39b) and (39d), we can
equivalently rewrite Eq. (39) as

(40a) $1-p_{0}+p_{1}-p_{2}>0$ (40b) $1-p_{0}>0,$ (40c) $1+p_{0}>0,$ (40d)
$1+p_{0}+p_{1}+p_{2}>0,$ (40e) $1-p_{1}+p_{0}p_{2}-p_{0}^{2}>0.$

Actually, under the condition of $\cos\theta\neq 1$, we can first obtain

(41a) $\displaystyle
1-p_{0}+p_{1}-p_{2}=(2-s_{1})(2-s_{2})(1+\cos\theta)+s_{1}s_{2}\omega_{0}(1-\cos\theta)>0,$
(41b) $1+p_{0}=1+(s_{1}-1)(1-s_{2})>0,$ (41c)
$1-p_{0}=1-(s_{1}-1)(1-s_{2})>0,$ (41d)
$1+p_{0}+p_{1}+p_{2}=s_{2}(1-\cos\theta)(2-s_{1})(1-\omega_{0})>0,$

where $0<\omega_{0}<1$ and $0<s_{1},s_{2}<2$ have been used.

To prove Eq. (40e), we first introduce the parameters $A$ and $B$,

(42) $\displaystyle A$
$\displaystyle=s_{1}(1-s_{2})(2-s_{1})+\omega_{0}s_{2}(1-s_{1})(2-s_{2}),$
$\displaystyle B$
$\displaystyle=s_{1}s_{2}(s_{1}+s_{2}-s_{1}s_{2})=s_{1}s_{2}[1-(1-s_{1})(1-s_{2})]>0,$

and can express the part on the left hand side of Eq. (40e) as

(43) $\displaystyle 1-p_{0}^{2}-p_{1}+p_{0}p_{2}$
$\displaystyle=[s_{1}(1-s_{2})(2-s_{1})+\omega_{0}s_{2}(1-s_{1})(2-s_{2})](1-\cos\theta)$
$\displaystyle+s_{1}s_{2}(s_{1}+s_{2}-s_{1}s_{2})=A(1-\cos\theta)+B.$

(i) If $0<s_{1}\leq 1$ and $0<s_{2}\leq 1$, we can obtain $A\geq 0$ and the
following equation,

(44) $\displaystyle 1-p_{0}^{2}-p_{1}+p_{0}p_{2}=A(1-\cos\theta)+B>0.$

(ii) If $0<s_{1}\leq 1$ and $1<s_{2}<2$, we have

(45) $A=s_{1}(1-s_{2})(2-s_{1})+\omega_{0}s_{2}(1-s_{1})(2-s_{2})\geq
s_{1}(1-s_{2})(2-s_{1}),$

which can be used to derive

(46) $\displaystyle 1-p_{0}^{2}-p_{1}+p_{0}p_{2}\geq
2s_{1}(1-s_{2})(2-s_{1})+B=s_{1}(2-s_{2})[1+(1-s_{1})(1-s_{2})]>0.$

(iii) If $1<s_{1}<2$ and $0<s_{2}\leq 1$, one can obtain

(47)
$A=s_{1}(1-s_{2})(2-s_{1})+\omega_{0}s_{2}(1-s_{1})(2-s_{2})>s_{2}(1-s_{1})(2-s_{2}),$

then we have

(48) $\displaystyle
1-p_{0}^{2}-p_{1}+p_{0}p_{2}>2s_{2}(1-s_{1})(2-s_{2})+B=s_{2}(2-s_{1})[1+(1-s_{1})(1-s_{2})]>0.$

(iv) If $1<s_{1}<2$ and $1<s_{2}<2$, one can derive

(49)
$A=s_{1}(1-s_{2})(2-s_{1})+\omega_{0}s_{2}(1-s_{1})(2-s_{2})>s_{1}(1-s_{2})(2-s_{1})+s_{2}(1-s_{1})(2-s_{2}).$

With the help of above equation and let
$C=s_{1}(1-s_{2})(2-s_{1})+s_{2}(1-s_{1})(2-s_{2})$, we can obtian

(50) $\displaystyle 1-p_{0}^{2}-p_{1}+p_{0}p_{2}$ $\displaystyle>2A+B>2C+B$
$\displaystyle=(2-s_{1})(2-s_{2})[1+(s_{1}-1)(1-s_{2})]>0.$

Based on above results (i)-(iv), one can find that Eq. (40e) indeed holds
under the conditions of $0<\omega_{0}<1$ and $0<s_{1},s_{2}<2$. Thus the roots
of characteristic polynomial $p(\lambda)$ are located in the field
$|\lambda|<1$ under the condition of $\cos\theta\neq 1$.

Now let us focus on the case of $\cos\theta=1$. We note that the roots of the
characteristic polynomial $p(\lambda)$ are continuous functions of
$\cos\theta$, and hence the roots of characteristic polynomial satisfy the
condition $|\lambda_{k}|\leq 1$ ($k$=1, 2 and 3). In addition, we would also
like to point out that for the special case of $\cos\theta=1$, one can adopt
the reductive approach [77] to obtain $|\lambda_{k}|\leq 1$ ($k$=1, 2 and 3).

From above discussion, we can find that the roots of characteristic polynomial
satisfy the condition $|\lambda_{k}|\leq 1$ ($k$=1, 2 and 3), thus the present
MRT-LB model based macroscopic numerical scheme (22) is unconditionally
stable.

##  5 Numerical results and discussion

To test the capacity of the MRT-LB model based macroscopic four-level
numerical scheme (22), we considered the diffusion equation (1) with the
following initial and boundary conditions [78],

(51) $\displaystyle\phi(x,0)=\sin(\pi x),~{}0\leq x\leq 1,$
$\displaystyle\phi(0,t)=\phi(1,t)=0,~{}t>0,$

and obtained the analytical solution of this problem,

(52) $\phi(x,t)=\sin(\pi x)e^{-\kappa\pi^{2}t}.$

In the implementation of the macroscopic four-level numerical scheme (22), the
analytical solution (52) is used to initialize the variable $\phi$ at first
three time levels. We first performed some simulations under different values
of the mesh Fourier number $\epsilon$ (or different diffusion coefficients for
the specified time and space steps), and plotted the results in Fig. 2 where
$\Delta x=0.025$, $\Delta t=0.01875$, the weight coefficient $\omega_{0}$ and
the relaxation parameters $s_{1}$ and $s_{2}$ are determined from Eqs. (30)
and (31). As shown in this figure, the numerical results are in good agreement
with the analytical solutions at time $t=12$.

Figure 2: The numerical and analytical solutions under different values of
discretization parameter $\epsilon$.

In addition, to quantitatively measure the deviation between the numerical
results and analytical solutions, the following root-mean-square error
($RMSE$) is adopted [17],

(53) $RMSE=\sqrt{\frac{\sum_{j=1}^{N}\big{[}\phi(j\Delta x,n\Delta
t)-\phi^{*}(j\Delta x,n\Delta t)\big{]}^{2}}{N}},$

where $N$ is the number of grid points, $\phi$ and $\phi^{*}$ are the
numerical and analytical solutions. Based on the definition of $RMSE$, one can
evaluate the convergence rate ($C_{R}$) of numerical scheme with the following
formula,

(54) $C_{R}=\log\big{(}RMSE_{\Delta x}/RMSE_{\Delta x/2}\big{)}/\log 2.$

We now focus on the $RMSE$ and convergence rate of the macroscopic four-level
numerical scheme (22). For this purpose, we conducted some simulations under
different values of space step $\Delta x$ and mesh Fourier number $\epsilon$
[see Table 1 for details, the weight coefficient $\omega_{0}$ and relaxation
parameters $s_{1}$ and $s_{2}$ are given by Eqs. (30) and (31)], and measured
the $RMSE$ between the analytical and numerical solutions at time $t=12$. As
seen from Table 2 and Fig. 3, the MRT-LB model based macroscopic four-level
numerical scheme indeed has a sixth-order convergence rate in space once the
conditions of Eqs. (30) and (31) are satisfied, and the numerical results with
a smaller mesh Fourier number $\epsilon$ are more accurate. However, as
demonstrated in the previous accuracy analysis, if the conditions of Eqs. (30)
and (31) are not met, the MRT-LB model based macroscopic numerical scheme
could not achieve the sixth-order accuracy in space. To confirm this
statement, we also carried out some simulations under the conditions of Eqs.
(30) and (31a) [Eq. (31b) is not satisfied], and presented the errors at
different space steps in Table 3 and Fig. 4. From these table and figure, one
can observe that the MRT-LB model based macroscopic four-level numerical
scheme is just fourth-order accurate in space. Moreover, if both Eqs. (31a)
and (31b) are not satisfied in our simulations, the MRT-LB model based
macroscopic four-level numerical scheme, as the commonly used LB method [28],
only has a second-order convergence rate, as reported in Table 4 and Fig. 5.
We note that these results are in agreement with our theoretical analysis.

Table 1: The values of some parameters in the MRT-LB based macroscopic four-level numerical scheme under the conditions of Eqs. (30) and (31). $\epsilon$ | $\omega_{0}$ | $s_{1}$ | $s_{2}$
---|---|---|---
$0.1$ | $0.8310204592587027$ | $0.9159290534201945$ | $1.1450386147380731$
$0.15$ | $0.8101626131270389$ | $0.775103705680168$ | $1.1476236168426883$
$0.175$ | $0.8370678725639358$ | $0.6352970255557769$ | $1.1776696173022918$
$0.2$ | $0.870066309422671$ | $0.49037716562528605$ | $1.2047312964902426$
$0.24$ | $0.9274277013170459$ | $0.2626707812024917$ | $1.2388413217086902$
Table 2: The $RMSE$ and $C_{R}$ of the MRT-LB based macroscopic four-level numerical scheme under the conditions of Eqs. (30) and (31). $\epsilon$ | $RMSE_{\Delta x=0.1}$ | $RMSE_{\Delta x=0.05}$ | $RMSE_{\Delta x=0.025}$ | $C_{R}$
---|---|---|---|---
$0.1$ | $8.59\times 10^{-10}$ | $1.42\times 10^{-11}$ | $2.57\times 10^{-13}$ | $\sim 5.85$
$0.15$ | $3.99\times 10^{-8}$ | $6.56\times 10^{-10}$ | $1.04\times 10^{-11}$ | $\sim 5.95$
$0.175$ | $1.19\times 10^{-7}$ | $1.95\times 10^{-9}$ | $3.11\times 10^{-11}$ | $\sim 5.95$
$0.2$ | $3.04\times 10^{-7}$ | $5.00\times 10^{-9}$ | $7.96\times 10^{-11}$ | $\sim 5.95$
$0.24$ | $1.31\times 10^{-6}$ | $2.15\times 10^{-8}$ | $3.43\times 10^{-10}$ | $\sim 5.95$
Figure 3: The convergence rate of the MRT-LB model based macroscopic four-level FD scheme under the conditions of Eqs. (30) and (31). Table 3: The $RMSE$ and $C_{R}$ of the MRT-LB based macroscopic four-level numerical scheme under the conditions of Eqs. (30) and (31a). $\epsilon$ | $\omega_{0}$ | $s_{1}$ | $s_{2}$ | $RMSE_{\Delta x=0.1}$ | $RMSE_{\Delta x=0.05}$ | $RMSE_{\Delta x=0.025}$ | $C_{R}$
---|---|---|---|---|---|---|---
$0.1$ | $0.8$ | $1$ | $12/11$ | $4.68\times 10^{-7}$ | $3.08\times 10^{-8}$ | $1.96\times 10^{-9}$ | $\sim 3.95$
$0.15$ | $0.7$ | $1$ | $42/41$ | $2.21\times 10^{-6}$ | $1.46\times 10^{-7}$ | $9.30\times 10^{-9}$ | $\sim 3.95$
$0.175$ | $0.65$ | $1$ | $78/79$ | $5.13\times 10^{-6}$ | $3.39\times 10^{-7}$ | $2.16\times 10^{-8}$ | $\sim 3.95$
$0.2$ | $0.6$ | $1$ | $18/19$ | $9.84\times 10^{-6}$ | $6.49\times 10^{-7}$ | $4.14\times 10^{-8}$ | $\sim 3.95$
$0.24$ | $0.52$ | $1$ | $78/89$ | $2.19\times 10^{-5}$ | $1.44\times 10^{-6}$ | $9.16\times 10^{-8}$ | $\sim 3.95$
Figure 4: The convergence rate of the MRT-LB model based macroscopic four-level FD scheme under the conditions of Eqs. (30) and (31a). Table 4: The $RMSE$ and $C_{R}$ of the MRT-LB based macroscopic four-level numerical scheme under the condition of Eq. (30) [Eqs. (31a) and (31b) are not satisfied]. $\epsilon$ | $\omega_{0}$ | $s_{1}$ | $s_{2}$ | $RMSE_{\Delta x=0.1}$ | $RMSE_{\Delta x=0.05}$ | $RMSE_{\Delta x=0.025}$ | $C_{R}$
---|---|---|---|---|---|---|---
$0.1$ | $0.8$ | $1$ | $1$ | $5.65\times 10^{-4}$ | $1.49\times 10^{-4}$ | $3.81\times 10^{-5}$ | $\sim 1.95$
$0.15$ | $0.7$ | $1$ | $1$ | $1.77\times 10^{-4}$ | $4.62\times 10^{-5}$ | $1.17\times 10^{-5}$ | $\sim 1.96$
$0.175$ | $0.65$ | $1$ | $1$ | $8.77\times 10^{-5}$ | $2.40\times 10^{-5}$ | $6.18\times 10^{-6}$ | $\sim 1.91$
$0.2$ | $0.6$ | $1$ | $1$ | $3.76\times 10^{-4}$ | $1.00\times 10^{-4}$ | $2.57\times 10^{-5}$ | $\sim 1.93$
$0.24$ | $0.52$ | $1$ | $1$ | $8.55\times 10^{-4}$ | $2.27\times 10^{-4}$ | $5.79\times 10^{-5}$ | $\sim 1.94$
Figure 5: The convergence rate of the MRT-LB model based macroscopic four-
level finite-difference scheme under the condition of (30).

##  6 Conclusions

In this work, we first developed a MRT-LB model for 1D diffusion equation
where the D1Q3 lattice structure is considered, and then obtained a mesoscopic
MRT-LB model based macroscopic four-level FD scheme. Through the theoretical
analysis, one can find that the macroscopic four-level numerical scheme is
unconditionally stable, and can achieve the sixth-order accuracy in space once
the weight coefficient $\omega_{0}$ and the relaxation parameters $s_{1}$ and
$s_{2}$ satisfy Eqs. (30) and (31). And also, if only the conditions of Eqs.
(30) and (31a) are met, the macroscopic four-level numerical scheme would have
a fourth-order convergence rate in space. Moreover, if only Eq. (30) is
satisfied, the macroscopic four-level numerical scheme would be second-order
accurate in space, which is the same as the commonly used LB method. In
addition, compared to the previous work [17], the present macroscopic
numerical scheme could be more accurate through adjusting the weight
coefficient $\omega_{0}$ and the relaxation parameters $s_{1}$ and $s_{2}$
properly.

We conducted some simulations to test the MRT-LB model based macroscopic four-
level FD scheme, and found that the numerical results agree well with our
theoretical analysis. Additionally, it is also interesting that the mesoscopic
two-level LB method can be equivalent to a macroscopic four-level FD scheme,
this indicates that the distribution function based mesoscopic method is more
efficient than the corresponding partial-differential-equation based
macroscopic numerical scheme, this may be because more information is included
in the distribution function.

Finally, we would like to point out that in the D1Q3 LB method based
macroscopic numerical schemes, the present macroscopic numerical scheme (22)
with the sixth-order convergence rate is the most accurate for 1D diffusion
equations. We also note that the present results are only limited to 1D case,
the MRT-LB model based macroscopic numerical schemes for two and three-
dimensional diffusion problems would be considered in a future work.

## References

  * [1] R. Benzi, S. Succi, M. Vergassola, The lattice Boltzmann equation: Theory and applications, Phys. Reports 222 (1992) 145-197.
  * [2] S. Chen, G. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329-364.
  * [3] C. K. Aidun J. R. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (2010) 439-472.
  * [4] L. Chen, Q. Kang, Y. Mu, Y.-L. He, W.-Q. Tao, A critical review of the pseudopotential multiphase lattice Boltzmann model: Methods and applications, Int. J. Heat Mass Transf. 76 (2014) 210-236.
  * [5] Q. Li, K. H. Luo, Q. J. Kang, Y. L. He, Q. Chen, Q. Liu, Lattice Boltzmann methods for multiphase flow and phase-change heat transfer, Prog. Energy Combus. Sci. 52 (2016) 62-105.
  * [6] H. Liu, Q. Kang, C. R. Leonardi, S. Schmieschek, A. Narváez, B. D. Jones, J. R. Williams, A. J. Valocchi, J. Harting, Multiphase lattice Boltzmann simulations for porous media applications, Comput. Geosci. 20 (2016) 777-805.
  * [7] H. Wang, X. Yuan, H. Liang, Z. Chai, B. Shi, A brief review of the phase-field-based lattice Boltzmann method for multiphase flows, Capillary 2 (2019) 33-52.
  * [8] K. V. Sharma, R. Straka, F. W. Tavares, Current status of lattice Boltzmann methods applied to aerodynamic, aeroacoustic, and thermal flows, Prog. Aero. Sci. 115 (2020) 100616.
  * [9] D. A. Wolf-Gladrow, Lattice-gas cellular automata and lattice Boltzmann models: An introduction, Springer, Berlin, 2000.
  * [10] S. Succi, The lattice Boltzmann equation for fluid dynamics and beyond, Oxford University Press, Oxford, 2001.
  * [11] M. C. Sukop, D. T. Thorne, Jr, Lattice Boltzmann modeling: An introduction for geoscientists and engineers, Springer, New York, 2006.
  * [12] Z. Guo, C. Shu, Lattice Boltzmann method and its applications in engineering, World Scientific Publishing Co. Pte. Ltd., Singapore, 2013.
  * [13] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E. M. Viggen, The lattice Boltzmann method: Principles and practice, Springer, Switzerland, 2017.
  * [14] S. P. Dawson, S. Chen, G. D. Doolen, Lattice Boltzmann computations for reaction-diffusion equations, J. Chem. Phys. 98 (1993) 1514-1523.
  * [15] D. A. Wolf-Gladrow, A lattice Boltzmann equation for diffusion, J. Stat. Phys. 79 (1995) 1023-1032.
  * [16] C. Huber, B. Chopard, M. Manga, A lattice Boltzmann model for coupled diffusion, J. Comput. Phys. 229 (2010) 7956-7976.
  * [17] S. Suga, An accurate multi-level finite difference scheme for 1D diffusion equations derived from the lattice Boltzmann method, J. Stat. Phys. 140 (2010) 494-503.
  * [18] P. V. Leemput, C. Vandekerckhove, W. Vanroose, D. Roose, Accuracy of hybrid lattie Boltzmann/finite difference schemes for reaction-diffusion systems, Multiscale Model. Simul. 6 (2007) 838-857.
  * [19] R. G. M. van der Sman, M. H. Ernst, Convection-diffusion lattice Boltzmann scheme for irregular lattices, J. Comput. Phys. 160 (2000) 766-782.
  * [20] I. Ginzburg, Equilibrium-type and link-type lattice Boltzmann models for generic advection and anisotropic-dispersion equation, Adv. Water Resour. 28 (2005) 1171-1195.
  * [21] I. Rasin, S. Succi, W. Miller, A multi-relaxation lattice kinetic method for passive scalar diffusion, J. Comput. Phys. 206 (2005) 453-462.
  * [22] I. Ginzburg, Multiple anisotropic collisions for advection-diffusion lattice Boltzmann schemes, Adv. Water Resour. 51 (2013) 381-404.
  * [23] B. Shi, Z. Guo, Lattice Boltzmann model for nonlinear convection-diffusion equations, Phys. Rev. E 79 (2009) 016701.
  * [24] B. Chopard, J. L. Falcone, J. Latt, The lattice Boltzmann advection-diffusion model revisited, Eur. Phys. J. Special Topics 171 (2009) 245-249.
  * [25] H. Yoshida, M. Nagaoka, Multiple-relaxation-time lattice Boltzmann model for the convection and anisotropic diffusion equation, J. Comput. Phys. 229 (2010) 7774-7795.
  * [26] Z. Chai, T. S. Zhao, Lattice Boltzmann model for the convection-diffusion equation, 87 (2013) 063309.
  * [27] R. Huang, H. Wu, A modified multiple-relaxation-time lattice Boltzmann model for convection-diffusion equation, J. Comput. Phys. 274 (2014) 50-63.
  * [28] Z. Chai, B. Shi, Z. Guo, A multiple-relaxation-time lattice Boltzmann model for general nonlinear anisotropic convection-diffusion equations, J. Sci. Comput. 69 (2016) 355-390.
  * [29] O. Aursjø, E. Jettestuen, J. L. Vinningland, A. Hiorth, An improved lattice Boltzmann method for simulating advective-diffusive processes in fluids, J. Comput. Phys. 332 (2017) 363-375.
  * [30] L. Li, Multiple-time-scaling lattice Boltzmann method for the convection diffusion equation, Phys. Rev. E 99 (2019) 063301.
  * [31] M. Hirabayashi, Y. Chen, H. Ohashi, The lattice BGK model for the Poisson equation, JSME Int. J. Ser. B 44 (2001) 45-52.
  * [32] Z. Chai, B. Shi, A novel lattice Boltzmann model for the Poisson equation, Appl. Math. Model. 32 (2008) 2050-2058.
  * [33] Z. Chai, H. Liang, R. Du, B. Shi, A lattice Boltzmann model for two-phase flow in porous media, SIAM J. Sci. Comput. 41 (2019) B746-B772.
  * [34] H. Lai, C. Ma, Lattice Boltzmann method for the generalized Kuramoto-Sivashinsky equation, Physica A 388 (2009) 1405-1412.
  * [35] H. Otomo, B. M. Boghosian, F. Dubois, Two complementary lattice-Boltzmann-based analyses for nonlinear systems, Physica A 486 (2017) 1000-1011.
  * [36] Z. Chai, N. He, Z. Guo, B. Shi, Lattice Boltzmann model for high-order nonlinear partial differential equations, Phys. Rev. E 97 (2018) 013304.
  * [37] L. Zhong, S. Feng, P. Dong, S. Gao, Lattice Boltzmann schemes for the nonlinear Schrödinger equation, Phys. Rev. E 74 (2006) 036704.
  * [38] S. Succi, Lattice Boltzmann method for quantum field theory, J. Phys. A 40 (2007) F559-F567.
  * [39] S. Palpacelli, S. Succi, Numerical validation of the quantum lattice Boltzmann scheme in two and three dimensions, Phys. Rev. E 75 (2007) 066704.
  * [40] B. Shi, Z. Guo, Lattice Boltzmann model for the one-dimensional nonlinear Dirac equation, Phys. Rev. E 79 (2009) 066704.
  * [41] Y. H. Qian, D. d’Humières, P. Lallemand, Lattice BGK models for Navier-Stokes equation, Europhys. Lett. 17 (1992) 479-484.
  * [42] D. d’Humières, Generalized lattice-Boltzmann equations, in: B.D. Shizgal, D.P. Weave (Eds.), Rarefied Gas Dynamics: Theory and Simulations, in: Prog. Astronaut. Aeronaut., Vol. 159, AIAA, Washington, DC, 1992, pp. 450-458.
  * [43] Z. Chai, B. Shi, Multiple-relaxation-time lattice Boltzmann method for the Navier-Stokes and nonlinear convection-diffusion equations: Modeling, analysis, and elements, Phys. Rev. E 102 (2020) 023306.
  * [44] P. Lallemand, L.-S. Luo, Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability, Phys. Rev. E 61 (2000) 6546-6562.
  * [45] C. Pan, L.-S. Luo, C. T. Miller, An evaluation of lattice Boltzmann schemes for porous medium flow simulation, Comput. Fluids 35 (2006) 898-909.
  * [46] S. Cui, N. Hong, B. Shi, Z. Chai, Discrete effect on the halfway bounce-back boundary condition of multiple-relaxation-time lattice Boltzmann model for convection-diffusion equations, Phys. Rev. E 93 (2016) 043311.
  * [47] L.-S. Luo, W. Liao, X. Chen, Y. Peng, W. Zhang, Numerics of the lattice Boltzmann method: Effects of collision models on the lattice Boltzmann simulations, Phys. Rev. E 83 (2011) 056701.
  * [48] S. Chapman, T. G. Cowling, The mathematical theory of non-uniform gases, Cambridge University Press, Cambridge, 1970.
  * [49] E. Ikenberry, C. Truesdell, On the Pressures and the Flux of Energy in a Gas according to Maxwell’s Kinetic Theory, I, J. Rational Mech. Anal. 5 (1956) 1-54.
  * [50] W.-A. Yong, W. Zhao, L.-S. Luo, Theory of the lattice Boltzmann method: Derivation of macroscopic equations via the Maxwell iteration, Phys. Rev. E 93 (2016) 033310.
  * [51] D. J. Holdych, D. R. Noble, J. G. Georgiadis, R. O. Buckius, Truncation error analysis of lattice Boltzmann methods, J. Comput. Phys. 193 (2004) 595-619.
  * [52] A. Wagner, Thermodynamic consistency of liquid-gas lattice Boltzmann simulations, Phys. Rev. E 74 (2006) 056703\.
  * [53] F. Dubois, Equivalent partial differential equations of a lattice Boltzmann scheme, Comput. Math. Appl. 55 (2008) 1441-1449.
  * [54] F. Dubois, Third order equivalent equation of lattice Boltzmann scheme, Discret. Contin. Dyn. Syst. 23 (2009) 221-248.
  * [55] F. Dubois, P. Lallemand, Towards higher order lattice Boltzmann schemes,J. Stat. Mech. 2009 (2009) P06006.
  * [56] D. d’Humières, I. Ginzburg, Viscosity independent numerical errors for lattice Boltzmann models: From recurrence equations to ”magic” collision numbers, Comput. Math. Appl. 58 (2009) 823-840.
  * [57] I. Ginzburg, Truncation errors, exact and heuristic stability analysis of two-relaxation-times lattice Boltzmann schemes for anisotropic advection-diffusion equation, Commun. Comput. Phys. 11 (2012) 1439- 1502\.
  * [58] M. Junk, A. Klar, and L.-S. Luo, Asymptotic analysis of the lattice Boltzmann equation, J. Comput. Phys. 201 (2005) 676-704.
  * [59] M. Junk, A finite difference interpretation of the lattice Boltzmann method, Numer. Meth. Part. Diff. Equ. 17 (2001) 383-402.
  * [60] T. Inamuro, A lattice kinetic scheme for incompressible viscous flows with heat transfer, Phil. Trans. R. Soc. Lond. A 360 (2002) 477-484.
  * [61] M. G. Ancona, Fully-Lagrangian and lattice-Boltzmann methods for solving systems of conservation equations, J. Comput. Phys. 115 (1994) 107-120.
  * [62] E. C. Du Fort and S. P. Frankel, Stability conditions in the numerical treatment of parabolic differential equations, Math. Comput. 7 (1953) 135-152.
  * [63] X. He, Q. Zou, L.-S. Luo, M. Dembo, Analytic solutions of simple flows and analysis of nonslip bundary conditions for the lattice Boltzmann BGK model, J. Stat. Phys. 87 (1997) 115-136.
  * [64] Z. Guo, C. Zheng, Analysis of lattice Boltzmann equation for microscale gas flows: Relaxation times, boundary conditions and the Knudsen layer, Int. J. Comput. Fluid Dyn. 22 (2008) 465-473.
  * [65] Z. Chai, B. Shi, Z. Guo, J. Lu, Gas flow through square array of circular cylinders with Klinkerberg effect: A lattice Boltzmann study, Commun. Comput. Phys. 8 (2010) 1052-1073.
  * [66] Q. Li, Z. Zheng, S. Wang, and J. Liu, A multilevel finite difference scheme for one-dimensional Burgers equation derived from the lattice Boltzmann method, J. Appl. Math. 2012 (2012) 925920.
  * [67] Z. Chai, C. Huang, B. Shi, Z. Guo, A comparative study on the lattice Boltzmann models for predicting effective diffusivity of porous media, Int. J. Heat Mass Transf. 98 (2016) 687-696.
  * [68] J. Latt, B. Chopard, Lattice Boltzmann method with regularized pre-collision distribution functions, Math. Comput. Simul. 72 (2016) 165-168
  * [69] L. Wang, B. Shi, Z. Chai, Regularized lattice Boltzmann model for a class of convection-diffusion equations, Phys. Rev. E 92 (2015) 043311.
  * [70] X. Yang, B. Shi, Z. Chai, Generalized modification in the lattice Bhatnagar-Gross-Krook model for incompressible Navier-Stokes equations and convection-diffusion equations, Phys. Rev. E 90 (2014) 013309.
  * [71] L. Wang, J. Mi, X. Meng, Z. Guo, A localized mass-conserving lattice Boltzmann approach for non-Newtonian fluid flows, Commun. Comput. Phys. 17 (2015) 908-924.
  * [72] L. Wang, W. Zhao, X.-D. Wang, Lattice kinetic scheme for the Navier-Stokes equations coupled with convection-diffusion equations, Phys. Rev. E 98 (2018) 033308.
  * [73] J. W. Thomas, Numerical Partial Differential Equations: Finite Difference Methods, Springer, New York, 1995.
  * [74] E. J. Routh, A treatise of the stability of a given state of motion, Macmillan, London, 1877.
  * [75] A. Hurwitz, Ueber die Bedingungen, unter welchen eine Gleichung nur Wurzeln mit negativen reellen Theilen besitzt, Math. Ann. 1895, 46: 273-284.
  * [76] F. R. Gantmacher, Applications of the theory of matrices, Interscience Publishers, Inc., New York, 1959.
  * [77] J. H. Miller, On the location of zeros of certain classes of polynomials with applications to numerical analysis, J. Inst. Maths Applics 1971, 8: 397-406.
  * [78] C. A. J. Fletcher, Computational techniques for fluid dynamics 1: Fundamental and general techniques, second edition, Springer, Berlin, 1991.


# A Human-Centered Approach for
Bootstrapping Causal Graph Creation

Minh Q. Tram<EMAIL_ADDRESS>The University of Texas at
ArlingtonArlingtonTXUSA76019 ,  Nolan B. Gutierrez
<EMAIL_ADDRESS>The University of Texas at
ArlingtonArlingtonTXUSA76019 and  William J. Beksi<EMAIL_ADDRESS>The
University of Texas at ArlingtonArlingtonTXUSA76019

###### Abstract.

Causal inference, a cornerstone in disciplines such as economics, genomics,
and medicine, is increasingly being recognized as fundamental to advancing the
field of robotics. In particular, the ability to reason about cause and effect
from observational data is crucial for robust generalization in robotic
systems. However, the construction of a causal graphical model, a mechanism
for representing causal relations, presents an immense challenge. Currently, a
nuanced grasp of causal inference, coupled with an understanding of causal
relationships, must be manually programmed into a causal graphical model. To
address this difficulty, we present initial results towards a human-centered
augmented reality framework for creating causal graphical models. Concretely,
our system bootstraps the causal discovery process by involving humans in
selecting variables, establishing relationships, performing interventions,
generating counterfactual explanations, and evaluating the resulting causal
graph at every step. We highlight the potential of our framework via a
physical robot manipulator on a pick-and-place task.

Human-Robot Interaction, Augmented Reality, Causal Graphs

##  1\. Introduction

Causality allows for deciphering complex relationships and drawing informed
conclusions. This foundational understanding is beginning to make significant
contributions in the field of robotics, reshaping the path of its advancement
(Li, 2023). However, constructing accurate causal models, including causal
inference and graph representations, is a daunting task. This stems from the
inherent intricacy of the interactions, often compounded by a myriad of
variables, which often remain elusive or hidden in the environment.

A causal graphical model, aka causal graph, serves as a visual mathematical
representation of the relationships between variables in a system. Delineating
the direction and structure of these relationships allows for identifying
potential causes and effects, and disentangling confounding variables
(Schölkopf et al., 2021). By employing causal graphs, robots can predict the
outcomes of their actions more accurately, make informed decisions in dynamic
environments, and adapt more fluidly to new situations with a deeper
comprehension of the underlying mechanisms (Xiong et al., 2016; Lee et al.,
2021; Cao et al., 2021).

Parallel to these advancements, the use of simulations and virtual interfaces
offers a more economical means of gathering diverse and valuable data, which
is especially conducive to improving machine learning algorithms (Mittal et
al., 2023). For example, the integration of virtual, augmented, and mixed
reality (VAMR) has shown promising results in fostering a more intuitive and
enhanced interaction between human operators and robots (e.g., (Puljiz et al.,
2019, 2020; Tram et al., 2023)). VAMR interfaces provide a way for operators
to have a more immersive control experience, bridging the gap between the
digital and physical realms.

Fig. 1. The point of view of the operator from the AR overlaid physical
workspace. The operator can interact with the robot and provide additional
context to aid the robot’s understanding via on-the-fly construction of causal
graphs.

In this work we present preliminary results on a new framework that leverages
VAMR technologies to address the challenge of causal graph construction, Fig.
1. To summarize, our contributions are the following.

  * •

We create an augmented reality (AR) interface that allows an operator to
naturally construct a causal graph on the fly.

  * •

We provide a human-centered approach to highlight critical information when
establishing a causal graph of a scene.

  * •

Our system ensures that a robot prioritizes human insights during the initial
stages of scene comprehension.

##  2\. Background and Motivation

###  2.1. Causal Graphs and Causal Inference

VBT Fig. 2. A causal graph representing the relationship between a robot’s
battery level $B$, the terrain roughness $T$, and the robot’s velocity $V$
during navigation. The solid edge represents a well-established cause and
effect, while the dashed edge represents an indirect or latent confounding
variable.

To establish a foundation for causal inference, Pearl introduced a number of
mathematical formulations, notably structural causal models (SCMs),
intervention and do-calculus, and counterfactuals, to rigorously define the
relationships between variables. For example, consider a mobile robot
navigation scenario. The causal graph that represents the relationship between
the robot’s velocity, the battery level, and the roughness or difficulty of
the terrain can be represented as a directed acyclic graph (DAG), Fig. 2. In
this scenario, we can represent the relationships between the variables using
an SCM $f$ as

(1) $V=f(B,T,U_{V}),$

where $V$ is the velocity, $B$ is the battery level, $T$ is the terrain
roughness, and $U_{V}$ is an internal variable that is unobserved, but has an
effect on the velocity of the robot. An intervention in this situation would
be to charge the robot’s battery regardless of any other factor hence ensuring
it is always full. This can be represented using do-calculus, i.e.,

(2) $P(V\,|\,do(B=\text{full}))=\sum_{T}P(V\,|\,B=\text{full},T=t)P(T=t).$

In (2), estimates of the speed of the robot are given through the intervention
by averaging over all possible values of $T$. Assuming that there are no
confounding relationships between the variables, we have effectively broken
the natural causal relationship between $B$ and $V$, while leaving the causal
relationship between $T$ and $V$ intact. This is where counterfactuals can
help infer the causal relationship between $B$ and $V$ during an intervention,
which can aid in reinforcing or disproving causal relationships between
variables. This is important for both pre-planning, on-the-fly decision
making, and insights toward post-incident analysis. By using this mathematical
framework to understand the causal relationships between variables, we can
make informed predictions about a system. For instance, in the robot
navigation scenario we can use these relationships to pinpoint the root cause
of performance issues, pre-charge the robot’s battery if it is expected to go
through rough terrain that will drain the battery, or adjust the robot’s
behavior to achieve a desired outcome.

###  2.2. Motivation

The use of causal inference is not new in robotics (Brawer et al., 2020; Ho et
al., 2020; Diehl and Ramirez-Amaro, 2023; Ding et al., 2023). Nonetheless,
automating the creation of a viable causal graph remains an underexplored area
of research. Constructing a causal graph is hard due to the following reasons.

  * •

Complexity of the environment: Variables and their causal relationships must
be identified in dynamic environments while accounting for confounding
variables.

  * •

Causation vs. correlation: The distinction between causation and correlation
from observational data is not always clear, and definitive proof of causation
often requires more data.

  * •

High dimensionality: The number of variables in an environment can be large,
which can lead to complex and/or multiple causal graphs.

##  3\. Human-Centered Causal Graphical Models

In this section, we delineate the specific methodologies and procedures
employed to actualize human-centered causal graphs. First, we present a brief
overview of our framework’s components. Then, we provide a detailed
description of each component and the communication pipeline between them.

###  3.1. Component Overview and Consideration

Fig. 3. The visualization and interaction pipeline of our proposed framework.
The operator can interact with the robot using an overlay interface via
various modes and provide context hinting directly onto the workspace to aid
the robot in its understanding of the scene.

Our framework to facilitate robotic understanding through human-centered
causal graph construction is composed of three primary components: (i) an AR-
capable headset, (ii) a software interface to visualize the robot’s perception
and allow the operator to interact with the robot, and (iii) a control
interface to realize the operator’s commands and interventions. Each of these
components plays a pivotal role in the seamless execution of the system,
providing both unique and complementary functionalities. Fig. 3 details the
overall system architecture and the transmission pipeline between the
components.

For the AR headset, we opted for the Microsoft HoloLens 2 due to its
unparalleled capabilities. Its commercial availability ensures ease of access
and its extensive API support, via Microsoft’s Mixed Reality Toolkit (MRTK)
(mrt, 2024), offers a versatile platform for development and integration. The
HoloLens 2 is a self-contained head-mounted display (HMD) that is capable of
spatial mapping and tracking, while also providing a projected holographic
display. It permits the operator to operate in a hands-free, untethered
manner. This feature is crucial for our work for two primary reasons. First,
being able to superimpose additional information directly onto the workspace
enhances the operator’s understanding of the robot’s perception of the
environment. This is vital when visualizing the implications of the causal
graph on the robot’s physical actions in real-time. Second, the ability to
operate in a hands-free manner through different modalities (i.e., direct-hand
input, eye-tracking, and voice commands) allows the operator to swiftly and
naturally interact with both the robot’s actions and the underlying causal
structure. This ensures that the operator can manipulate the causal graph,
adjust parameters, and interact with robot perception and environment
understanding without the constraints of traditional input methods.

###  3.2. Interaction and Visualization Modalities

Relying on the HoloLens 2’s hand-tracking capabilities, our system allows the
operator to interact with the robot and environment through a number of
mechanisms. This includes direct manipulation of the robot’s joint positions
via virtual sliders and command buttons, providing a planning context such as
desired end-effector positions and approach vectors, and the ability to query
forward and inverse kinematics solutions. Furthermore, compared to traditional
VR devices where operator interactions are bound to handheld controllers and
the operator’s vision is obstructed by a display, the HoloLens 2 offers
greater interaction and visualization capabilities.

###  3.3. Causal Graph Construction and Intervention

In our human-centered framework, the focus is primarily on the dynamic
construction and intervention of causal graphs. To this end, we leveraged
Unity (uni, 2024) to render causal structures as 3D objects in the AR space
allowing the operator to interact with the causal graph directly. These
graphical constructs, composed of nodes and edges, offer a visual guide that
aids the operator in deciphering the complex interplay between robot actions
and environmental cues.

Nodes, visualized as tangible entities within the AR landscape, symbolize
specific events or actions that the robot can detect or perform. Their design
incorporates distinct visual features, such as shades, opacities, and
pulsations, communicating their present state or relevance. On the other hand,
edges provide visual insights into the direction, strength, or likelihood of
causal linkages. These are often characterized by variations in their
thickness, texture, or color gradient, allowing users to quickly comprehend
causal dependencies.

Within the immersive world facilitated by the HMD, operators can actively
modify and intervene in the causal relationships by adding, removing, or
modifying the nodes and edges of the current causal structure. The device’s
advanced hand-tracking feature supports intuitive gestures, such as pinching,
swiping, or rotating, enabling users to adjust causality strength, node
prominence, or even the direction of a causal link. Furthermore, the mobility
of the HMD lets an operator freely move around the environment and provide
context hints directly onto the robot workspace, accommodating seamless
adjustments and experiments with the causal graph.

###  3.4. Physical Robot Actuation and Simulation

The physical robot actuation and simulation components of our system, driven
by ROS 2 Humble (ros, 2024), facilitate communication between the robot and
Unity software layers. This allows the operator to interact with the robot and
environment through the AR interface. Concretely, it enables the translation
of intricate directives from a causal graph into precise physical movements of
the robot, ensuring adaptability and dexterity in diverse real-world
scenarios.

##  4\. Evaluation

###  4.1. Causal Graph Creation

Using a UR5e robot arm provisioned with a Robotiq 2F-85 gripper and tasked
with picking up a finite set of known objects and placing them at a desired
location, we construct a causal graph to represent the relationship between
the robot’s actions and the environment. To achieve this, the potential
variables related to our task are represented as follows.

  1. (1)

$Robot$: The robot state.

  2. (2)

$Gripper$: The gripper state.

  3. (3)

$Target$: The robot’s motion planning to achieve the goal position.

  4. (4)

$Goal$: The desired object placement location.

  5. (5)

$Type$: The object’s type.

  6. (6)

$Weight$: The object’s weight (light - heavy).

  7. (7)

$Size$: The object’s size (small - large).

  8. (8)

$Texture$: The object’s texture (smooth - rough).

  9. (9)

$Rigidity$: The object’s rigidness (soft - hard).

  10. (10)

$Success$: The success rate of a pick-and-place sequence.

For the purpose of demonstration, we construct this causal graph under the
assumption that external factors (e.g., robot power loss, environmental
disruptions, etc.) do not occur and are insulated from the causal
representation. As shown in Fig. 4, the relationships between the variables
can be visually represented as a graph for an easier understanding of
potential causality and the interplay between nodes. Within the graph, each
edge represents a direct causal relationship between the robot and the
operating environment. For example, if the weight of the object changes then
it could influence the behavior of the gripper (e.g., the gripper force or
grasp vector).

Our framework can directly create the causal graph in the AR space, Fig. 5.
The operator interacts via the AR interface to label the nodes and edges, set
the content of each node, and adjust the strength of each edge. The operator
can also communicate with the robot and the environment through the AR
interface to obtain additional information and context to aid in the
construction of the causal graph. The causal graph can then be parsed and
interpreted by the motion planning backend during the planning step to ensure
that the robot’s actions are consistent with the state of the graph. Algorithm
1 demonstrates how the robot consumes a causal graph generated by the operator
during each motion planning step.

WeightSizeTextureRigidityTypeGripperTargetRobotGoalSuccess Fig. 4. An example
of the relationships between the operator-identified variables for a pick-and-
place task. Each edge corresponds to a direct causal relationship between the
connected variables. Changes in the source of an edge have a direct
consequence on the target of an edge. Algorithm 1 Human-in-the-loop Causal
Graph Iteration

1:Human has initialized the causal graph

2:Causal graph is DAG, all nodes have properties

3:while true do

4: if operator is modifying causal graph then

5: Wait

6: else

7: Execute with guidance from causal graph

8: Check execution status

9: if status is SUCCESS then

10: Continue

11: else

12: Notify Operator

13: Ask Continue?

14: end if

15: end if

16:end while

Fig. 5. The operator interacting with the robot and environment through the AR
interface to construct a causal graph. Nodes and edges are visualized and can
be manipulated directly in the AR space.

###  4.2. Causal Graph Intervention

For a pick-and-place task, we use our constructed causal graph to evaluate
whether the relationships between the defined variables are valid, which in
turn can be inferred by the robot’s performance. For instance, if the robot
consistently fails to pick up a cordless drill but it has no issues with other
objects, then we can be relatively certain that several known properties or
relationships (e.g., weight, texture, etc.) may be detrimentally affecting the
robot. With this human-in-the-loop approach, we can quickly devise an
intervention procedure to decide if certain variables are causing degraded
performance. For example, we can intervene by setting the object’s weight to a
specific value and then benchmark the robot’s performance under the assumption
that the robot did not account for the object’s weight properly.

If the robot’s performance is influenced by the weight of the object, then we
can make adjustments to the robot’s behavior to ensure optimal performance
(e.g., adjusting the pick points, etc.). Alternatively, should the
intervention be inconclusive, then we can intervene by fixing other object
properties to a specific value and then benchmark the robot’s performance.
However, further investigation may be required to understand the exact reason.
For example, the gripper might not be suited for the cordless drill’s shape,
or the perception system may not be identifying the drill’s correct pick
points. Since our framework has the capability to adjust causal graph
structure on the fly, we can quickly update the graph to reflect the new
structure thus allowing the robot to make more informed decisions.

###  4.3. Discussion

Our framework provides advantages over traditional methods for constructing
causal graphs. The AR interface allows an operator to interact with the robot
and the environment directly and more intuitively, making it easier for humans
to understand the system’s dynamics. It also presents a visual understanding
of causal relationships, which allows operators to make informed decisions on
how to optimize the system. For example, the operator can visualize in real
time where interventions might be needed, or where potential failures could
occur.

Constructing causal graphs this way introduces a level of modularity and
scalability. For instance, causal elements can be easily adjusted to fit the
needs of the system or the entire graph can be replaced by constructing a
replacement from scratch. Yet, this methodology is not without its
limitations. A criticism of causal graphs is that they can be an
oversimplification of the environment and interplay between variables. Real-
world robotic operations may involve numerous subtle interactions, not all of
which can be accurately captured in a graph.

##  5\. Conclusion

This paper provided preliminary results on a human-centered approach for
automating the construction of causal graphical models. Our framework combines
the strengths of VAMR and simulation technologies to address the dilemma of
bootstrapping the creation of a causal graph. Furthermore, we highlighted the
ability to visualize, intervene, and update causal graphs on a physical robot
for a real-world pick-and-place task. In the future we will work on addressing
the challenges of oversimplification and iteratively updating causal graphs.

## References

  * (1)
  * mrt (2024) 2024\.  _Microsoft Mixed Reality Toolkit_.  https://github.com/microsoft/MixedRealityToolkit-Unity
  * ros (2024) 2024\.  _Robot Operating System_.  https://www.ros.org
  * uni (2024) 2024\.  _Unity 3D Game Engine_.  https://unity.com
  * Brawer et al. (2020) Jake Brawer, Meiying Qin, and Brian Scassellati. 2020.  A Causal Approach to Tool Affordance Learning. In _Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems_. 8394–8399.
  * Cao et al. (2021) Yu Cao, Boyang Li, Qian Li, Adam Stokes, David Ingram, and Aristides Kiprakis. 2021\.  Reasoning Operational Decisions for Robots via Time Series Causal Inference. In _Proceedings of the IEEE International Conference on Robotics and Automation_. 6124–6131.
  * Diehl and Ramirez-Amaro (2023) Maximilian Diehl and Karinne Ramirez-Amaro. 2023.  A Causal-Based Approach to Explain, Predict and Prevent Failures in Robotic Tasks.  _Robotics and Autonomous Systems_ 162 (2023), 104376.
  * Ding et al. (2023) Wenhao Ding, Haohong Lin, Bo Li, and Ding Zhao. 2023\.  CausalAF: Causal Autoregressive Flow for Safety-Critical Driving Scenario Generation. In _Proceedings of the Conference on Robot Learning_. PMLR, 812–823.
  * Ho et al. (2020) Seng-Beng Ho, Mark Edmonds, and Song-Chun Zhu. 2020\.  Actional-Perceptual Causality: Concepts and Inductive Learning for AI and Robotics. In _Proceedings of the IEEE Symposium Series on Computational Intelligence_. 442–448.
  * Lee et al. (2021) Tabitha E. Lee, Jialiang Alan Zhao, Amrita S. Sawhney, Siddharth Girdhar, and Oliver Kroemer. 2021\.  Causal Reasoning in Simulation for Structure and Transfer Learning of Robot Manipulation Policies. In _Proceedings of the IEEE International Conference on Robotics and Automation_. 4776–4782.
  * Li (2023) Yangming Li. 2023\.  Deep Causal Learning for Robotic Intelligence.  _Frontiers in Neurorobotics_ 17 (2023).
  * Mittal et al. (2023) Mayank Mittal, Calvin Yu, Qinxi Yu, Jingzhou Liu, Nikita Rudin, David Hoeller, Jia Lin Yuan, Ritvik Singh, Yunrong Guo, Hammad Mazhar, Ajay Mandlekar, Buck Babich, Gavriel State, Marco Hutter, and Animesh Garg. 2023.  Orbit: A Unified Simulation Framework for Interactive Robot Learning Environments.  _IEEE Robotics and Automation Letters_ (2023).
  * Puljiz et al. (2020) David Puljiz, Franziska Krebs, Fabian Bosing, and Bjorn Hein. 2020\.  What the HoloLens Maps Is Your Workspace: Fast Mapping and Set-up of Robot Cells via Head Mounted Displays and Augmented Reality. In _Proceeding of the IEEE/RSJ International Conference on Intelligent Robots and Systems_. 11445–11451.
  * Puljiz et al. (2019) David Puljiz, Erik Stöhr, Katharina S. Riesterer, Björn Hein, and Torsten Kröger. 2019.  General Hand Guidance Framework Using Microsoft HoloLens. In _Proceeding of the IEEE/RSJ International Conference on Intelligent Robots and Systems_. 5185–5190.
  * Schölkopf et al. (2021) Bernhard Schölkopf, Francesco Locatello, Stefan Bauer, Nan Rosemary Ke, Nal Kalchbrenner, Anirudh Goyal, and Yoshua Bengio. 2021\.  Toward Causal Representation Learning.  _Proc. IEEE_ 109, 5 (2021), 612–634.
  * Tram et al. (2023) Minh Q Tram, Joseph M Cloud, and William J Beksi. 2023\.  Intuitive Robot Integration via Virtual Reality Workspaces. In _Proceedings of the IEEE International Conference on Robotics and Automation_. 11654–11660.
  * Xiong et al. (2016) Caiming Xiong, Nishant Shukla, Wenlong Xiong, and Song-Chun Zhu. 2016.  Robot Learning with a Spatial, Temporal, and Causal and-or Graph. In _Proceedings of the IEEE International Conference on Robotics and Automation_.