#  A Dominant Interferer-based Approximation for Uplink SINR Meta Distribution
in Cellular Networks

Yujie Qin, Mustafa A. Kishk, Member, IEEE, and Mohamed-Slim Alouini, Fellow,
IEEE Yujie Qin and Mohamed-Slim Alouini are with Computer, Electrical and
Mathematical Sciences and Engineering (CEMSE) Division, King Abdullah
University of Science and Technology (KAUST), Thuwal, 23955-6900, Saudi
Arabia. Mustafa Kishk is with the Department of Electronic Engineering,
Maynooth University, Maynooth, W23 F2H6, Ireland. (e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>slim.alouini@kaust.edu.sa).

###### Abstract

This work studies the signal-to-interference-plus-noise-ratio (SINR) meta
distribution for the uplink transmission of a Poisson network with Rayleigh
fading by using the dominant interferer-based approximation. The proposed
approach relies on computing the mix of exact and mean-field analysis of
interference. In particular, it requires the distance distribution of the
nearest interferer and the conditional average of the rest of the
interference. Using the widely studied fractional path-loss inversion power
control and modeling the spatial locations of base stations (BSs) by a Poisson
point process (PPP), we obtain the meta distribution based on the proposed
method and compare it with the traditional beta approximation, as well as the
exact results obtained via Monte-Carlo simulations. Our numerical results
validate that the proposed method shows good matching and is time competitive.

###### Index Terms:

Meta distribution, approximation, stochastic geometry, Poisson network,
reliability, uplink

##  I Introduction

For a wireless network, the system’s performance depends heavily on the
spatial configuration of nodes. The accurate modeling of the locations of
nodes provides key system insights on the network performance [1]. Stochastic
geometry provides tractable mathematical tools to capture the randomness of
the nodes [2]. While investigating the spatial average performance metrics,
such as coverage probability, is important, the performance of an individual
link can no longer be ignored, as it reveals the reliability and quality of
service of a network [3]. On the side of operators, obtaining information
about the percentile of users achieving a certain performance measure is a
fundamental design objective. This high-level system insight is so-called
SINR/SIR meta distribution, defined as a complementary cumulative distribution
function (CCDF) of the conditional success probability (probability that the
user equipment (UE) achieves cellular coverage, conditioned on the realization
of point processes), and used to investigate the discrepancies among UE [4].

The exact expression of SINR/SIR meta distribution is typically computed by
using Gil-Pelaez theorem [5], however, it is generally hard to obtain as it
requires imaginary moments. Some inequalities, such as Markov, Chebyshev, and
Chernoff, were provided in [3] to bound the meta distribution. A comprehensive
introduction about meta distribution, as well as some examples on PPP networks
and Poisson bipolar networks, was provided in [6, 7]. Specifically, in [7],
the author provided a closed-form result of SIR meta distribution for PPP and
Poisson bipolar network by using nearest-interferer-only approximation.
Authors in [8] computed the SIR meta distribution for any independent fading
and analyzed the separability. Authors in [9] approximated the meta
distribution of a general network by shifting the meta distribution of a
Poisson network.

For investigating SIR/SINR meta distribution in different scenarios, authors
in [10] studied the Poisson line Cox bipolar networks and Poisson stick line
Cox bipolar networks; authors in [11] analyzed cellular networks with BS
cooperation; authors in [12] characterized a $K$-tier HetNet combined by PPP
and MCP; and authors in [13] considered a typical cell. Besides, the SINR meta
distribution of UAV-involved networks was studied in [14, 15].

In our previous work [16], we proposed a dominant interferer-based
approximation as an alternative method of computing downlink meta distribution
in Poisson networks with Rayleigh fading, and it shows good matching with the
exact value. In this work, we extend it into uplink scenarios and complete the
analysis of the proposed approximation. Different from downlink scenarios, the
system models of the uplink case are more complex, for instance, the locations
of interferers are modeled by a non-homogeneous PPP, the nearest interferer
may be closer than the typical UE to the BS, and the transmit power is a
function of the distance to the serving BS.

##  II System Model

This work focuses on the uplink transmission of a PPP wireless cellular
network. The locations of BSs are distributed according to a PPP, denoted by
$\Phi$ with density $\lambda$. Assuming BSs split the resources to serve the
UE, therefore, the UE point process is generated by independently and
uniformly selecting one UE in each Voronoi cell formed by $\Phi$. From the
perspective of a BS, the locations of interferers can be approximately modeled
as a non-homogeneous PPP, $\Phi_{i}$, with intensity function
$\lambda_{i}(d)=\lambda(1-\exp(-\pi\lambda d^{2}))$, where $d$ denotes the
distance to the BS [17]. Here, we condition a BS to be at the origin and this
BS becomes the typical BS when averaging over the point process. We focus on
the SINR of the typical BS which is equivalent to any other arbitrary
deterministic location owing to the stationarity of PPP.

When UE is served by its nearest BS, a fractional path-loss inversion power
control with compensation factor $\epsilon$ and a standard path-loss model
with exponent $\alpha>2$ are used. Besides, a Rayleigh fading channel model is
used in this work.

By using the inversion power control technology [18, 4], the truncated
transmit power of a user in the uplink transmission is given by

$\displaystyle p_{t}(R_{u})=\left\\{\begin{aligned} &\rho
R_{u}^{\alpha\epsilon},\quad\text{if}\quad
0<R_{u}<\bigg{(}\frac{p_{\max}}{\rho}\bigg{)}^{\frac{1}{\alpha\epsilon}},\\\
&p_{\rm max},\quad\text{if}\quad
R_{u}\geq\bigg{(}\frac{p_{\max}}{\rho}\bigg{)}^{\frac{1}{\alpha\epsilon}},\\\
\end{aligned}\right.$ (1)

where $R_{u}$ is the Euclidean distance between the user and the BS, $\rho$ is
a power control parameter to adjust the received power at the serving BS,
$\epsilon$ is the compensation factor, and $p_{\rm max}$ is the maximum
transmit power. Consequently, the received power at the BS is

$\displaystyle p_{r}=Hp_{t}(R_{u})R_{u}^{-\alpha},$ (2)

in which $H$ present the fading gain which follows exponential distribution
with mean of unity.

The uplink SINR of the typical BS is

$\displaystyle{\rm
SINR}=\frac{Hp_{r}}{\sigma^{2}+I}=\frac{Hp_{t}(R_{u})R_{u}^{-\alpha}}{\sigma^{2}+I},$
(3)

where $\sigma^{2}$ is the thermal noise and $I$ is the aggregate interference,

$\displaystyle I$
$\displaystyle=\sum_{x\in\Phi_{i}}H_{x}p_{t,x}(R_{i,x})D_{i,x}^{-\alpha},$ (4)

in which $R_{i,x}$ is the distance between the interferer $x$ to its serving
BS, $H_{x}$, $p_{r,x}$ and $D_{i,x}$ denote the channel fading, transmit power
and distances of the interferer $x$, respectively. In what follows, the
conditional success probability of the typical BS is given by

$\displaystyle P_{s}(\theta)=\mathbb{P}({\rm SINR}>\theta\mid\Phi_{i}).$ (5)

For an arbitrary realization of $\Phi_{i}$, we analyze the SINR meta
distribution, which is defined as the fraction of links achieving the uplink
coverage (SINR above $\theta$), and its argument denotes the system
reliability.

######  Definition 1 (SINR Meta Distribution).

The SINR meta distribution of the uplink is defined in [3] as

$\displaystyle\bar{F}_{P_{s}}(\theta,\gamma)=\mathbb{P}(P_{s}(\theta)>\gamma),$
(6)

where $\gamma\in[0,1]$.

##  III Mathematical Analysis

In this section, we provide the analysis for the dominant interfer-based
approximation for SINR meta distribution in the case of uplink transmission.
To compare the proposed approximation with the traditional method, we also
derive the beta approximation, which is based on the first and the second
moment of the conditional success probability. To do so, we first derive some
important distance distributions.

Recall that $R_{u}$ denotes the distance to the serving BS, $R_{i}$ and
$D_{i}$ are the distances between the interferer to its serving BS and typical
BS, respectively, and $D_{1}$ is the distance from the typical BS to the first
nearest UE $D_{1}=\min(D_{i,x})$, e.g., $D_{1}=D_{i,x_{1}}$.

######  Lemma 1 (Distance Distribution).

The probability density function (PDF) and cumulative distribution function
(CDF) of $R_{u}$ is given by

$\displaystyle f_{R_{u}}(r)$ $\displaystyle=2\pi\lambda\exp(-\lambda\pi
r^{2}),\quad r\geq 0,$ $\displaystyle F_{R_{u}}(r)$
$\displaystyle=1-\exp(-\lambda\pi r^{2}),\quad r\geq 0.$ (7)

Based on the strongest association police, $R_{i}$ is upper bounded by
$D_{i}$, which forms a truncated Rayleigh distribution, and given by

$\displaystyle f_{R_{i}}(r|D_{i})=\frac{2\pi\lambda\exp(-\lambda\pi
r^{2})}{1-\exp(-\lambda\pi D_{i}^{2})},\quad 0\leq r\leq D_{i}.$ (8)

The PDF of $D_{1}$ is given by

$\displaystyle f_{D_{1}}(r)=$ $\displaystyle 2\pi\lambda(1-\exp(-\pi\lambda
r^{2}))r$ $\displaystyle\times\exp(-2\pi\lambda\int_{0}^{r}(1-\exp(-\pi\lambda
z^{2}))z{\rm d}z).$ (9)

In what follows, we compute the SINR meta distribution based on the proposed
dominant interferer-based approximation and beta approximation, which are
provided in the next two subsections.

###  III-A Proposed Approximation

As mentioned in [16], the closer interferer has higher impact on the system
performance and dominant interferer-based approximation works well in downlink
transmission. Similarly, we follow the same method and rewrite the
interference term by considering the nearest interfering UE exactly while the
rest of the interferers in average. To simplify the notation, we use $R_{i,1}$
and $p_{t,1}(R_{i,1})$ to represent $R_{i,x_{1}}$ and
$p_{t,x_{1}}(R_{i,x_{1}})$, and $p_{t}^{\prime}(r)=\min(p_{\rm max},\rho
r^{-\alpha\epsilon})$.

######  Lemma 2 (Approximated Interference).

The approximation of the interference is

$\displaystyle I_{1}\approx Hp_{t,1}(R_{i,1})R_{1}^{-\alpha}+G(D_{1}),$ (10)

in which

$\displaystyle G(D_{1})=$ $\displaystyle
2\pi\lambda\int_{D_{1}}^{\infty}\int_{0}^{z}p_{t}^{\prime}(r)2\pi\lambda\exp(-\pi\lambda
r^{2})z^{-\alpha+1}{\rm d}r{\rm d}z.$ (11)

###### Proof:

See Appendix -A. ∎

Based on the approximated aggregate interference, the approximation of the
conditional success probability is given in the following lemma.

######  Lemma 3 (Approximated Conditional Success Probability).

The conditional success probability is approximated by

$\displaystyle\mathbb{P}({\rm SINR}>\theta\mid R_{u})$
$\displaystyle\approx\exp\bigg{(}\theta(\sigma^{2}+G(D_{1}))R_{u}^{\prime-1}\bigg{)}\frac{1}{1+\theta
p_{t}(D_{1})D_{1}^{-\alpha}R_{u}^{\prime-1}},$ (12)

in which $R_{u}^{\prime}=p_{t}(R_{u})R_{u}^{-\alpha}$.

###### Proof:

See Appendix -B. ∎

Now we are able to proceed to the final expression of the proposed
approximation in uplink scenarios.

######  Theorem 1 (Approximated SINR Meta Distribution).

The SINR meta distribution is approximated by

$\displaystyle\bar{F}_{P_{s}}^{\prime}(\theta,\gamma)\approx\int_{0}^{\infty}F_{R_{u}^{\prime}}(K(r,\theta,\gamma))f_{D_{1}}(r){\rm
d}r,$ (13)

in which

$\displaystyle K(r,\theta,\gamma)$
$\displaystyle\approx-\frac{1}{\kappa(r)}+\frac{1}{S(r)}W\bigg{(}0,\frac{S(r)\exp(S(r)\kappa^{-1}(r))}{\gamma\kappa(r)}\bigg{)},$
(14)

where $\kappa(r)=\theta\bar{p}_{t}(r)r^{-\alpha}$ and
$S(r)=\theta(\sigma^{2}+G(r))$,

$\displaystyle\bar{p}_{t}(r)$
$\displaystyle=\int_{0}^{r}p_{t}^{\prime}(r)f_{R_{i}}(x|r){\rm d}x,$
$\displaystyle F_{R_{u}^{\prime}}(x)$ $\displaystyle=\left\\{\begin{aligned}
&F_{R_{u}}((x\rho)^{\frac{1}{(1-\epsilon)\alpha}}),\quad\text{if}\quad
x<p_{\rm max}\bigg{(}\frac{p_{\rm
max}}{\rho}\bigg{)}^{\frac{1}{\epsilon\alpha}-\alpha},\\\ &F_{R_{u}}((xp_{\rm
max})^{\frac{1}{\alpha}}),\quad\text{if}\quad x>p_{\rm
max}\bigg{(}\frac{p_{\rm
max}}{\rho}\bigg{)}^{\frac{1}{\epsilon\alpha}-\alpha},\\\
\end{aligned}\right.$

where $F_{R_{u}}(r)$ is the CDF of $R_{u}$ given in (7) Lemma 1.

###### Proof:

See Appendix -C. ∎

###  III-B Beta Approximation

To compare the proposed approximation with existing method, in this
subsection, we derive the expression of beta approximation, which requires the
first and the second moment of the network.

######  Lemma 4 (Laplace Transform of the Interference).

The Laplace transform of the interference is given by

$\displaystyle\mathcal{L}_{I}(s)$
$\displaystyle=\exp\bigg{(}-2\pi\lambda\int_{0}^{\infty}\int_{0}^{r}\frac{2\pi\lambda\exp(-\lambda\pi
r^{2})}{1+(sp_{t}^{\prime}(x))^{-1}r^{\alpha}}{\rm d}xr{\rm d}r\bigg{)}.$ (15)

###### Proof:

See Appendix -D. ∎

Next, we compute the $b$-th moment of the network.

######  Lemma 5 ($b$-th Moment).

The $b$-th moment of the conditional success probability is

$\displaystyle M_{b}(\theta)$
$\displaystyle=\int_{0}^{\infty}f_{R_{u}}(z)\exp\bigg{(}-\kappa_{b}(z)\sigma^{2}\bigg{)}\mathcal{L}_{I}(\kappa_{b}(z)){\rm
d}z,$ (16)

where $\kappa_{b}(r)=b\frac{\theta r^{\alpha}}{p_{t}(r)}$.

###### Proof:

See Appendix -E. ∎

Typically, (6) is solved by using Gil-Pelaez theorem [5], which requires the
imaginary moments,

$\displaystyle\bar{F}_{P_{s}}(\theta,\gamma)=\frac{1}{2}+\frac{1}{\pi}\int_{0}^{\infty}\frac{1}{t}\Im(\exp(-jt\log\gamma)M_{jt}(\theta)){\rm
d}t,$ (17)

where $j$ is the imaginary unit and $M_{jt}(\theta)$ is obtained by replacing
the $b$ in the $b$-th moment by $jt$ and $\Im(\cdot)$ is the imaginary part of
a complex number.

Generally, computing (17) is tricky since it requires the imaginary moments.
Alternatively, we can use beta approximation to compute, which only requires
the first and the second moment and highly improves the analysis tractability,
and this approximation shows good matching for a wide range of scenarios [3].

######  Remark 1 (Beta Approximation).

The beta approximation is given by

$\displaystyle\bar{F}^{\prime\prime}_{P_{s}}(\theta,\gamma)\approx 1-$
$\displaystyle
I_{\gamma}\bigg{(}\frac{M_{1}(\theta)(M_{1}(\theta)-M_{2}(\theta))}{M_{2}(\theta)-M_{1}^{2}(\theta)},$
$\displaystyle\frac{(M_{1}(\theta)-M_{2}(\theta))(1-M_{1}(\theta))}{M_{2}(\theta)-M_{1}^{2}(\theta)}\bigg{)},$
(18)

where $I_{x}(a,b)=\frac{\int_{0}^{x}t^{a-1}(1-t)^{b-1}{\rm d}t}{B(a,b)}$, and
$B(a,b)=\int_{0}^{1}t^{a-1}(1-t)^{b-1}{\rm d}t$.

##  IV Numerical Results

(a) $\epsilon=0.4$ (b) $\epsilon=0.4$ (c) $\epsilon=0.8$ (d) $\epsilon=0.8$

Figure 1:  SINR meta distribution of uplink transmission of PPP networks. The
solid lines are exact values based on simulations, dash lines are for beta
approximation, and markers for the proposed approximation, respectively.

Figure 2:  SINR meta distribution of uplink transmission of PPP networks at
(a) $\gamma=0.6$ and (b) $\gamma=0.9$. The solid lines are exact values based
on simulations, dash lines are for beta approximation, and markers for the
proposed approximation, respectively.

In this section, we validate the proposed approximation and compare it with
the beta approximation and the exact value of the SINR meta distribution via
Monte-Carlo simulations with a large number of iterations to ensure accuracy.
For the simulation of the considered setup, we first generated two
realizations of PPPs for the locations of BSs and UE. UE is located in the
Voronoi cells formed by BSs and in each resource block only one UE is active.
Next, we compute the conditional success probability of each link (for each
fixed realization, channel fading varies for $10^{4}$ times to obtain the
conditional success probability) and obtain the CCDF of the conditional
success probability. Finally, we repeat this process for $10^{2}$ times.
Unless stated otherwise, we use the system parameters listed herein. The
maximum transmit power of UE is $p_{\rm max}=200$ mW, the power control factor
is $\rho=0.008$ mW, the noise power is set as $\sigma^{2}=10^{-9}$ W, and the
path-loss exponent is $\alpha=4$.

In Fig. 1, we plot the exact value, beta approximation, and the proposed
approximation of SINR meta distribution against $\theta$ or $\gamma$ at
different values of $\gamma$ or $\theta$ under $\epsilon=0.4,0.8$,
respectively. The proposed approximation shows good matching at all values of
$\gamma$ and $\theta$. Slightly different from the downlink scenario in [16],
the results of the proposed approximation is very close to beta approximation,
and difficult to tell which approximation has a better performance. For the
time of computing one point, the proposed approximation is slightly shorter
than the beta approximation. Even though the proposed approximation does not
provide a dramatic improvement in performance, we still believe it is
meaningful. While the beta approximation behaves more like a fitting of meta
distribution, as mentioned in [3], the reason for using beta distribution is
that its high order moments fit well with the moments of the network, and some
other approximations, such as nearest neighbor-only approximation, Markov and
Chebyshev bounds, show gaps with SINR meta distribution, the proposed method
provides a good matching approximation to the meta distribution. That is, in
some special cases, e.g., Rayleigh fading models, (6) can be solved by using
the dominant interferer-based approximation, and from this point of view, the
proposed approximation completes the analysis of meta distribution.

Fig. 2 provides some insights about the effect of $\epsilon$ on the proposed
approximation and SINR meta distribution. With the increases of $\epsilon$,
the system reliability first increases slowly, then increases dramatically
until the maximum achievable value which is limited by the maximum transmit
power and interference. This is because of the transmit power of UE and
distance to the serving BS. While the increase of $\epsilon$ increases the
transmit power of UE and compensates for the path-loss, it also increases the
interference. Therefore, an optimal $\epsilon$ exists to maximize the system’s
reliability. Besides, we notice that the proposed approximation, as well as
the beta approximation, shows a small gap with the exact value of meta
distribution at low values of $\theta$.

##  V Conclusion

This work analyzes a dominant interferer-based approximation in uplink
networks. Compared with the downlink scenario, the uplink analysis is more
complex since the truncated power inversion control is used and the locations
of interferers follow a non-homogeneous PPP. We first show that the proposed
approximation is operable in an uplink scenario which follows a similar way as
downlink scenarios. To validate the accuracy, we compare it with the
traditional method, beta approximation, as well as the exact value obtained
via simulation. The proposed approximation shows a good matching results in
uplink scenarios for different values of the path-loss compensation factor.
Consequently, this work completes the analysis of the dominant interferer-
based approximation of meta distribution and proves that the proposed
approximation works well in both uplink and downlink scenarios.

###  -A Proof of Lemma 2

Similar to [16, Lemma 2], we consider the nearest interferer exactly and the
rest interferers in mean sense. With that been said, conditioned on the
location of the nearest interferer, we compute the mean of the rest of
interferers. Consequently, the aggregate interference is approximated by the
exact term of nearest interferer and the conditional mean of the rest of the
terms,

$\displaystyle I_{1}$
$\displaystyle=\sum_{x\in\Phi_{i}}H_{x}p_{t,x}(R_{i,x})D_{i,x}^{-\alpha}$
$\displaystyle=H_{x_{1}}p_{t,x_{1}}(R_{i,x_{1}})D_{i,x_{1}}^{-\alpha}+\sum_{x\in\Phi_{i}\setminus
x_{1}}H_{x}p_{t,x}(R_{i,x})D_{i,x}^{-\alpha}$ $\displaystyle\approx
H_{x_{1}}p_{t,x_{1}}(R_{i,x_{1}})D_{i,x_{1}}^{-\alpha}+G(D_{i,x_{1}}),$ (19)

in which

$\displaystyle G(D_{i,x_{1}})$
$\displaystyle=\mathbb{E}\bigg{[}\sum_{x\in\Phi_{i}\setminus
x_{1}}H_{x}p_{t,x}(R_{i,x})D_{i,x}^{-\alpha}\bigg{]}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\mathbb{E}\bigg{[}\sum_{x\in\Phi_{i}\setminus
x_{1}}p_{t,x}(R_{i,x})D_{i,x}^{-\alpha}\bigg{]}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}2\pi\lambda\int_{D_{i,x_{1}}}^{\infty}\int_{0}^{z}(1-\exp(-\pi\lambda
z^{2}))$ $\displaystyle\quad\times\min(p_{\rm max},\rho
r^{\alpha\eta})f_{R_{u}}(r|z)(r)z^{-\alpha+1}{\rm d}r{\rm d}z,$ (20)

in which step (a) follows from the fact that all fading gains are i.i.d.
distributed with mean of unity and step (b) follows from Campbell’s theorem
[19] with conversion from Cartesian to polar coordinates.

###  -B Proof of Lemma 3

By using the approximated interference given in the previous lemma, the
conditional success probability is approximated by,

$\displaystyle\mathbb{P}({\rm SINR}>\theta\mid
R_{u})=\mathbb{P}\bigg{(}\frac{Hp_{t}(R_{u})R_{u}^{-\alpha}}{\sigma^{2}+I_{1}}>\theta\mid
R_{u}\bigg{)}$
$\displaystyle=\mathbb{P}\bigg{(}H>\frac{\theta(\sigma^{2}+I_{1})}{p_{t}(R_{u})R_{u}^{-\alpha}}\mid
R_{u}\bigg{)}$
$\displaystyle\approx\mathbb{P}\bigg{(}H>\theta(\sigma^{2}+H_{x_{1}}p_{t}(R_{i,1})D_{1}^{-\alpha}+G(D_{1}))\frac{R_{u}^{\alpha}}{p_{t}(R_{u})}\mid
R_{u}\bigg{)}$
$\displaystyle=\mathbb{E}\bigg{[}\exp\bigg{(}\theta(\sigma^{2}+H_{x_{1}}p_{t}(R_{i,1})D_{1}^{-\alpha}+G(D_{1}))\frac{R_{u}^{\alpha}}{p_{t}(R_{u})}\bigg{)}\bigg{]}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\exp\bigg{(}\theta(\sigma^{2}+G(D_{1}))R_{u}^{\prime-1}\bigg{)}\frac{1}{1+\theta
p_{t}(D_{1})D_{1}^{-\alpha}R_{u}^{\prime-1}},$ (21)

where step (a) follows from the use of MGF of exponential distribution.

###  -C Proof of Theorem 1

Recall that SINR meta distribution is the CCDF of the conditional success
probability, therefore, approximated by

$\displaystyle\bar{F}_{P_{s}}^{\prime}(\theta,\gamma)=\mathbb{P}(P_{s}(\theta)>\gamma)$
$\displaystyle\approx\mathbb{P}\bigg{(}\exp\bigg{(}\theta(\sigma^{2}+G(D_{1}))p_{t}^{-1}(R_{u})R_{u}^{\alpha}\bigg{)}$
$\displaystyle\quad\times\frac{1}{1+\theta
p_{t}(R_{i,1})D_{1}^{-\alpha}p_{t}^{-1}(R_{u})R_{u}^{\alpha}}>\gamma\bigg{)}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\mathbb{E}\bigg{[}\mathbb{P}\bigg{(}p_{t}(R_{u})R_{u}^{-\alpha}<-\frac{1}{\theta
p_{t}(R_{i,1})D_{1}^{-\alpha}}$
$\displaystyle\quad+\frac{1}{S(R_{1})}W\bigg{(}0,\frac{S(R_{i,1})\exp(S(D_{1})\theta^{-1}p_{t}^{-1}(R_{i,1})D_{1}^{\alpha})}{\gamma\theta
p_{t}(R_{i,1})R_{1}^{-\alpha}}\bigg{)}\bigg{)}\bigg{]}$
$\displaystyle=\int_{0}^{\infty}F_{R_{u}^{\prime}}(K(r,\theta,\gamma))f_{D_{1}}(r){\rm
d}r,$ (22)

in which step (a) can be solved by using simple algebraic manipulations and
the definition of the Lambert W function: $W(0,x)\exp(W(0,x))=x$, and

$\displaystyle K(r,\theta,\gamma)=-\frac{1}{\theta
p_{t}(R_{i,x_{1}})r^{-\alpha}}$
$\displaystyle\quad+\frac{1}{S(r)}W\bigg{(}0,\frac{S(r)\exp(S(r)\theta^{-1}p_{t}^{-1}(R_{i,1})r^{\alpha})}{\gamma\theta
p_{t}(R_{i,1})r^{-\alpha}}\bigg{)}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\approx}}-\frac{1}{\theta\bar{p}_{t}(r)r^{-\alpha}}+\frac{1}{S(r)}W\bigg{(}0,\frac{S(r)\exp(S(r)\theta^{-1}\bar{p}_{t}^{-1}(r)r^{\alpha})}{\gamma\theta\bar{p}_{t}(r)r^{-\alpha}}\bigg{)},$
(23)

in which step (a) is approximated by taking the integration of $p_{t}(r)$
separately.

Finally, the CDF of $R_{u}^{\prime}$ is computed by

$\displaystyle
F_{R_{u}^{\prime}}(x)=\mathbb{P}(R_{u}^{\prime}<x)=\mathbb{P}(p_{t}^{\prime}(R_{u})R_{u}^{-\alpha}<x)$
$\displaystyle=\left\\{\begin{aligned} &\mathbb{P}(\rho
R_{u}^{(1-\alpha)\epsilon}<x),\quad\text{if}\quad
0<R_{u}<\bigg{(}\frac{p_{\max}}{\rho}\bigg{)}^{\frac{1}{\alpha\epsilon}},\\\
&\mathbb{P}(p_{\rm max}R_{u}^{-\alpha}<x),\quad\text{if}\quad
R_{u}\geq\bigg{(}\frac{p_{\max}}{\rho}\bigg{)}^{\frac{1}{\alpha\epsilon}},\\\
\end{aligned}\right.$ (24)

proof completes by using $F_{R_{u}}(r)$.

###  -D Proof of Lemma 4

The Laplace transform of the aggregate interference is given by

$\displaystyle\mathcal{L}_{I}(s)=\mathbb{E}_{I}[\exp(-sI)]$
$\displaystyle=\mathbb{E}_{I}\bigg{[}\exp\bigg{(}-s\sum_{x\in\Phi_{i}}H_{x}p_{t,x}(R_{i,x})D_{i,x}^{-\alpha}\bigg{)}\bigg{]}$
$\displaystyle=\mathbb{E}_{I}\bigg{[}\prod_{x\in\Phi_{i}}\exp\bigg{(}-sH_{x}p_{t,x}(R_{i,x})D_{i,x}^{-\alpha}\bigg{)}\bigg{]}$
$\displaystyle=\mathbb{E}_{I}\bigg{[}\prod_{x\in\Phi_{i}}\mathbb{E}\bigg{[}\frac{1}{1+sp_{t,x}(R_{i,x})D_{i,x}^{-\alpha}}\bigg{]}\bigg{]}$
$\displaystyle=\mathbb{E}_{I}\bigg{[}\prod_{x\in\Phi_{i}}\int_{0}^{D_{i}}\frac{1}{1+sp_{t}^{\prime}(x)D_{i,x}^{-\alpha}}f_{R_{i}}(x|D_{i}){\rm
d}x\bigg{]}$
$\displaystyle=\exp\bigg{(}-2\pi\lambda\int_{0}^{\infty}(1-\exp(-\pi\lambda
r^{2}))$
$\displaystyle\times\bigg{(}1-\int_{0}^{r}\frac{1}{1+sp_{t}^{\prime}(x)r^{-\alpha}}f_{R_{i}}(x|r){\rm
d}x\bigg{)}r{\rm d}r\bigg{)}$
$\displaystyle=\exp\bigg{(}-2\pi\lambda\int_{0}^{\infty}\int_{0}^{r}\frac{2\pi\lambda\exp(-\lambda\pi
r^{2})}{1+(sp_{t}^{\prime}(x))^{-1}r^{\alpha}}{\rm d}xr{\rm d}r\bigg{)}$ (25)

###  -E Proof of Lemma 5

The $b$-th moment is computed by taking the expectation of the $b$-th moment
of the conditional success probability,

$\displaystyle M_{b}(\theta)$ $\displaystyle=\mathbb{E}[P_{s}^{b}(\theta)],$
(26)

in which the exact expression of the conditional success probability
(different from (Lemma 3)) is given by

$\displaystyle P_{s}(\theta)$ $\displaystyle=\mathbb{P}({\rm
SINR}>\theta|R_{u})=\mathbb{P}\bigg{(}\frac{Hp_{t}(R_{u})R_{u}^{-\alpha}}{\sigma^{2}+I}>\theta|R_{u}\bigg{)}$
$\displaystyle=\mathbb{P}\bigg{(}H>\frac{\theta
R_{u}^{\alpha}}{p_{t}(R_{u})}(\sigma^{2}+I)|R_{u}\bigg{)}$
$\displaystyle=\mathbb{E}_{I}\bigg{[}\exp\bigg{(}-\frac{\theta
R_{u}^{\alpha}}{p_{t}(R_{u})}(\sigma^{2}+I)\bigg{)}\bigg{]}$
$\displaystyle=\exp\bigg{(}-\frac{\theta
R_{u}^{\alpha}}{p_{t}(R_{u})}\sigma^{2}\bigg{)}\mathcal{L}_{I}\bigg{(}\frac{\theta
R_{u}^{\alpha}}{p_{t}(R_{u})}\bigg{)}.$ (27)

## References

  * [1] H. ElSawy, E. Hossain, and M. Haenggi, “Stochastic geometry for modeling, analysis, and design of multi-tier and cognitive cellular wireless networks: A survey,” _IEEE Communications Surveys & Tutorials_, vol. 15, no. 3, pp. 996–1019, 2013.
  * [2] H. ElSawy, A. Sultan-Salem, M.-S. Alouini, and M. Z. Win, “Modeling and analysis of cellular networks using stochastic geometry: A tutorial,” _IEEE Communications Surveys & Tutorials_, vol. 19, no. 1, pp. 167–203, 2016.
  * [3] M. Haenggi, “The meta distribution of the SIR in Poisson bipolar and cellular networks,” _IEEE Transactions on Wireless Communications_ , vol. 15, no. 4, pp. 2577–2589, 2015.
  * [4] H. ElSawy and M.-S. Alouini, “On the meta distribution of coverage probability in uplink cellular networks,” _IEEE Communications Letters_ , vol. 21, no. 7, pp. 1625–1628, 2017.
  * [5] J. Gil-Pelaez, “Note on the inversion theorem,” _Biometrika_ , vol. 38, no. 3-4, pp. 481–482, 1951.
  * [6] M. Haenggi, “Meta distributions—part 1: Definition and examples,” _IEEE Communications Letters_ , vol. 25, no. 7, pp. 2089–2093, 2021.
  * [7] ——, “Meta distributions—part 2: Properties and interpretations,” _IEEE Communications Letters_ , vol. 25, no. 7, pp. 2094–2098, 2021.
  * [8] K. Feng and M. Haenggi, “Separability, asymptotics, and applications of the SIR meta distribution in cellular networks,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 7, pp. 4806–4816, 2020.
  * [9] R. K. Ganti and M. Haenggi, “Asymptotics and approximation of the SIR distribution in general cellular networks,” _IEEE Transactions on Wireless Communications_ , vol. 15, no. 3, pp. 2130–2143, 2015.
  * [10] J. P. Jeyaraj, M. Haenggi, A. H. Sakr, and H. Lu, “The transdimensional Poisson process for vehicular network analysis,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 12, pp. 8023–8038, 2021.
  * [11] Q. Cui, X. Yu, Y. Wang, and M. Haenggi, “The SIR meta distribution in Poisson cellular networks with base station cooperation,” _IEEE Transactions on Communications_ , vol. 66, no. 3, pp. 1234–1249, 2017.
  * [12] C. Saha, M. Afshang, and H. S. Dhillon, “Meta distribution of downlink SIR in a Poisson cluster process-based HetNet model,” _IEEE Wireless Communications Letters_ , vol. 9, no. 12, pp. 2144–2148, 2020.
  * [13] P. D. Mankar, H. S. Dhillon, and M. Haenggi, “Meta distribution analysis of the downlink SIR for the typical cell in a Poisson cellular network,” in _2019 IEEE Global Communications Conference (GLOBECOM)_. IEEE, 2019, pp. 1–6.
  * [14] Y. Qin, M. A. Kishk, and M.-S. Alouini, “On the uplink SINR meta distribution of UAV-assisted wireless networks,” _IEEE Wireless Communications Letters_ , vol. 12, no. 4, pp. 684–688, 2023.
  * [15] ——, “On the downlink SINR meta distribution of UAV-assisted wireless networks,” _IEEE Transactions on Communications_ , 2023, to appear.
  * [16] ——, “A dominant interferer plus mean field-based approximation for SINR meta distribution in wireless networks,” _IEEE Transactions on Communications_ , vol. 71, no. 6, pp. 3663–3678, 2023.
  * [17] J. G. Andrews, A. K. Gupta, and H. S. Dhillon, “A primer on cellular network analysis using stochastic geometry,” _arXiv preprint arXiv:1604.03183_ , 2016\.
  * [18] H. ElSawy and E. Hossain, “On stochastic geometry modeling of cellular uplink transmission with truncated channel inversion power control,” _IEEE Transactions on Wireless Communications_ , vol. 13, no. 8, pp. 4454–4469, 2014\.
  * [19] M. Haenggi, _Stochastic geometry for wireless networks_. Cambridge University Press, 2012.


# High-order semi-Lagrangian kinetic scheme for compressible turbulence

Dominik Wilde  Department of Mechanical Engineering, University of Siegen,
Paul-Bonatz-Straße 9-11, 57076 Siegen-Weidenau, Germany  Institute of
Technology, Resource and Energy-efficient Engineering (TREE), Bonn-Rhein-Sieg
University of Applied Sciences, Grantham-Allee 20, 53757 Sankt Augustin,
Germany  Andreas Krämer  Department of Mathematics and Computer Science, Freie
Universität Berlin, Arnimallee 6, 14195 Berlin, Germany  Dirk Reith  Institute
of Technology, Resource and Energy-efficient Engineering (TREE), Bonn-Rhein-
Sieg University of Applied Sciences, Grantham-Allee 20, 53757 Sankt Augustin,
Germany  Fraunhofer Institute for Algorithms and Scientific Computing (SCAI),
Schloss Birlinghoven, 53754 Sankt Augustin, Germany  Holger Foysi  Department
of Mechanical Engineering, University of Siegen, Paul-Bonatz-Straße 9-11,
57076 Siegen-Weidenau, Germany

###### Abstract

Turbulent compressible flows are traditionally simulated using explicit time
integrators applied to discretized versions of the Navier-Stokes equations.
However, the associated Courant-Friedrichs-Lewy condition severely restricts
the maximum time step size. Exploiting the Lagrangian nature of the Boltzmann
equation’s material derivative, we now introduce a feasible three-dimensional
semi-Lagrangian lattice Boltzmann method (SLLBM), which circumvents this
restriction. While many lattice Boltzmann methods for compressible flows were
restricted to two dimensions due to the enormous number of discrete velocities
in three dimensions, the SLLBM uses only 45 discrete velocities. Based on
compressible Taylor-Green vortex simulations we show that the new method
accurately captures shocks or shocklets as well as turbulence in 3D without
utilizing additional filtering or stabilizing techniques, even when the time
step sizes are up to two orders of magnitude larger compared to simulations in
the literature. Our new method therefore enables researchers for the first
time to study compressible turbulent flows by a fully explicit scheme, whose
range of admissible time step sizes is only dictated by physics, while being
decoupled from the spatial discretization.

††preprint: APS/123-QED

##  I Introduction

One major challenge in fluid dynamics is the study of compressible turbulent
flows, involving intrinsic as well as variable density compressibility effects
[1, 2, 3, 4, 5, 6, 7]. Applications range from aviation [8] or astrophysics
[9] to the investigation of canonical flows like boundary layers [10], channel
flow [11, 12, 13], mixing layers [14, 15, 16, 17, 18], jets and aeroacoustics
[19, 20, 21] or shock-turbulence interaction [22], to only mention a few
considering the vast literature available. These flows feature both solenoidal
and dilatational structures, which constantly interact and possibly cause
shock waves [23, 5].

Numerical simulations have become an indispensable tool to understand their
physics, and many studies exploring compressible turbulent flows have been
conducted using high-order compact finite difference, optimized dispersion-
relation preserving schemes [24, 25, 26, 27, 28, 19, 29, 30] for the spatial
derivatives, often combined with low-dispersion-dissipation Runge-Kutta
schemes for time-integration [31, 19, 32]. Although these methods provide
accurate results, the time steps are generally small [33], because of the
methods’ Eulerian time derivatives, which describe how the variables of
interest pass through fixed locations in the field. Thus, the admissible time
step sizes are tightly linked to spatial resolution. This issue is for many
discretizations linked to the Courant-Friedrichs-Lewy (CFL) condition,

$c{\delta_{t}}/{\delta_{x}}<\mathrm{CFL}_{\mathrm{max}},$ (1)

using linear stability theory, relating a characteristic velocity $c$ to the
spatial and temporal discretization intervals $\delta_{x}$ and $\delta_{t}$,
respectively (see [34], for example). Though implicit time integration schemes
often provide larger stability domains, their application can be unfeasible
for transient problems due to their computational cost. Explicit time
integration schemes with scheme-specific $\mathrm{CFL_{max}}$, by contrast,
enforce small time steps $\delta_{t}$ for high flow velocities, typically
occurring in many high-speed compressible flows. Another obvious way to
circumvent the CFL condition in Eq. (1) is to incorporate Lagrangian time
derivatives, which track the motion of the variables of interest moving
through the domain.
In practice, Semi-Lagrangian (SL) schemes are used instead, which provide a
viable alternative to the discretization of Eulerian time derivatives. SL
schemes discretize the Lagrangian solution by tracking the trajectories back
in time. The prefix “semi” indicates that the trajectories’ end points usually
do not coincide with the simulation grid points, which requires application of
an appropriate interpolation scheme. SL methods were successfully incorporated
in algorithms solving the Navier-Stokes equations [35], although tracking of
the fluid trajectories was often found to be cumbersome, introducing
additional errors [36]. The major advantage when using SL schemes in kinetic
theory is that the trajectories are _linear_ , resulting in cancellation of
the tracking error. Consequently, SL schemes were both applied to the Vlasov
equation [37, 38, 39] and to the Bhatnagar-Gross-Krook (BGK)-Boltzmann
equation [40, 41]. Recently, we introduced the semi-Lagrangian lattice
Boltzmann method (SLLBM) [42, 43] for compressible flows [44], which solves
the lattice Boltzmann equation using a high-order SL streaming step.

In this article, we explore the capabilities of the SLLBM for three-
dimensional compressible flows. Furthermore, we demonstrate that the SLLBM
remains stable for time step sizes that exceed typical CFL constraints of
Eulerian solvers by orders of magnitude. To yield a lean scheme, the SLLBM is
combined with state-of-the-art cubature rules for the velocity discretization
[45, 46, 47]. This combination proves capable of modeling compressible
turbulence with time steps that are at least one order of magnitude larger
than in standard Eulerian methods and _decouple_ the spatial from the temporal
discretization.

###  I.1 Background

We start with a critical look at the more specialized lattice Boltzmann method
(LBM) [48]. Despite the successes of the standard LBM in the computation of
multiphase [49], particle-laden [50], thermal [51], or turbulent flows [52,
53], compressible LBM [54] were overlooked for a relatively long time, but
regained attraction during the last decade [55, 56, 57, 58, 59, 60, 61, 62,
63, 64, 65, 66, 67, 68]. Let us recall the force-free BGK-Boltzmann equation

$\frac{\partial f}{\partial t}+\bm{\xi}\cdot\nabla
f=-\frac{1}{\lambda}\left(f-f^{\mathrm{eq}}\right),$ (2)

with the continuous distribution function $f$, the equilibrium distribution
function $f^{\mathrm{eq}}$, the particle velocity $\bm{\xi}$, and the
relaxation time $\lambda$. To discretize Eq. (2), the original LBM is based on
three key principles. First, the equilibrium distribution function
$f^{\mathrm{eq}}$ is polynomially expanded into a series of Hermite
polynomials $\mathcal{H}^{(n)}$, with expansion coefficients being the
equilibrium moments $\mathbf{a}^{(n)}_{\mathrm{eq}}$ [69],

$f^{\mathrm{eq}}\approx\omega(\bm{\xi})\sum_{n=0}^{N}\frac{1}{n!}{\bm{a}}^{(n)}_{\mathrm{eq}}:{\bm{\mathcal{H}}}^{(n)},$
(3)

where $N$ is the expansion order and $\omega(\bm{\xi})$ the weight function.
Since $\bm{a}^{(n)}_{\mathrm{eq}}$ and $\bm{\mathcal{H}}^{(n)}$ are symmetric
tensors of rank $n$, the product involves contraction to all $D^{n}$ scalar
components, depending on dimension $D$. Second, a Gauß-Hermite quadrature is
applied to the unbounded velocity space of the Boltzmann equation, leading to
discrete particle velocity sets [69]. The moments are then found by the
quadrature

$\bm{a}^{(n)}=\int_{\mathbb{R}^{D}}\omega(\bm{\xi})\frac{f}{\omega(\bm{\xi})}\bm{\mathcal{H}}^{(n)}(\bm{\xi})d\bm{\xi}=\sum_{i=0}^{Q-1}f_{i}\bm{\mathcal{H}}^{(n)}(\bm{\xi}_{i}).$
(4)

with the weighted discrete distribution functions
$f_{i}=w_{i}f(\bm{\xi}_{i})/\omega(\bm{\xi}_{i}).$ The combination of $Q$
discrete particle velocities $\bm{\xi}_{i}$ and weights $w_{i}$, the velocity
set, is usually derived by the Gauß-product rule applied to a one-dimensional
Gauß-Hermite quadrature. Third, the discrete Boltzmann equation is integrated
along characteristics with an inherent Lagrangian discretization of the
Boltzmann equation’s material derivative to obtain a stable numerical scheme
and second-order temporal convergence [70].

Unfortunately, the LBM in its original formulation is mainly restricted to
Cartesian grids and velocity discretizations that match the regular lattices.
The customary “D2Q9” based on second-order expansion in Eq. (3) is plagued by
a cubic error being proportional to the Mach number. Consequently,
compressible simulations either demand correction terms that annihilate the
errors and restore Galilean invariance [71, 64] or higher-order
discretizations of Eqs. (3) and (4) which render abscissae that reside off-
lattice. Utilization of such velocity sets therefore requires an efficient
off-lattice Boltzmann solver. Previous Eulerian off-lattice Boltzmann schemes
[72, 73, 74, 75, 76], like finite difference or finite volume LBM, would be
suited in principle. However, they sacrifice the Lagrangian time integration
along characteristics. Moreover, their time step is severely restricted by a
CFL condition, Eq. (1), with respect to the fast discrete particle velocities.

In contrast, the SLLBM preserves all of the aforementioned key principles of
the LBM but it also decisively extends its capabilities. In previous works
[42, 43, 44] we have shown that a high-order interpolation increases the
spatial order of the method and nihilates mass losses. Also, we have
demonstrated the unconditional stability of the advection step, when
incorporating Gauß-Lobatto-Chebyshev nodes for the interpolation up to third
order, and that the stability is practically not affected even with fourth
order. The flexibility in terms of meshing and velocity sets encouraged us to
search for efficient quadrature rules solving the weight function. This
research led us to long-established cubature rules [45, 46, 47], i.e.
multivariate quadratures, which are often used in Kalman filters, e.g., in
[77].

##  II Methodology

###  II.1 Compressible semi-Lagrangian LBM

The compressible SLLBM uses the established lattice Boltzmann equation with
the BGK collision operator [42]

$h_{i}(\mathbf{x}\\!+\\!\delta_{t}\bm{\xi}_{i},t\\!+\\!\delta_{t})=h_{i}(\mathbf{x},t)-\frac{1}{\tau}\\!\left(h_{i}\\!\left(\mathbf{x},t\right)-h_{i}^{\mathrm{eq}}\\!\left(\mathbf{x},t\right)\right).$
(5)

Here, $h_{i}$ denotes either $f_{i}$ or the second distribution function
$g_{i}$ related to the variable heat capacity ratio $\gamma$. The shifted
dimensionless relaxation parameter
$\tau=\lambda/\delta_{t}+0.5=\mu/(Pc_{s}\delta_{t})+0.5$ depends on dynamic
viscosity $\mu$, lattice speed of sound $c_{s}$, and pressure $P=\rho T$ with
density $\rho$ and temperature $T$. The discrete equilibrium distribution
function $f_{i}^{\mathrm{eq}}$ is

$f_{i}^{\mathrm{eq},N}(\mathbf{x},t)=w_{i}\sum_{n=0}^{N}\frac{1}{n!}\bm{a}^{(n)}_{eq}(\mathbf{x},t):\bm{\mathcal{H}}_{i}^{(n)},$
(6)

and $g^{\mathrm{eq}}=(2C_{v}-D)Tf_{i}^{\mathrm{eq}}$, with heat capacity at
constant volume $C_{v}$, and number of dimensions $D$. Both
$\bm{a}^{(n)}_{eq}$ and
$\bm{\mathcal{H}}_{i}^{(n)}:=\bm{\mathcal{H}}^{(n)}(\bm{\xi}_{i})$ are listed
in Appendix A and B, respectively. To adjust the heat conductivity, a quasi-
equilibrium approach [78] is applied to Eq. (5), for more details see Appendix
D.

Density $\rho$, momentum $\rho u$, and energy $E$ are determined by

$\rho\\!=\\!\sum_{i=0}^{Q-1}f_{i},\;\rho
u\\!=\\!\sum_{i=0}^{Q-1}f_{i}\bm{\xi}_{i},\;2\rho
E=\sum_{i=0}^{Q-1}\left(f_{i}|\bm{\xi}_{i}|^{2}\\!+g_{i}\right).$ (7)

Appendix C demonstrates the connection of the present SLLBM model with the
macroscopic equations using a Chapman-Enskog analysis [79]. The integration
along characteristics, hidden in Eq. (5), incorporates a second-order temporal
error, whose order can be increased by multistep schemes [80]. In standard
LBMs, the particle velocities $\bm{\xi}_{i}$ in Eq. (5) are designed to end up
on one of the neighboring nodes, and the time step size is invariably set to
unity for the same purpose. By contrast, the SLLBM’s particle distribution
functions are still integrated along characteristics, but the departure points
may reside offside the grid, i.e. they are off-lattice. To recover the off-
lattice values an interpolation is needed. While several interpolation
strategies are possible, we chose a cell-oriented approach, which means that
once a departure point is identified, the degrees of freedom points
$\hat{h}_{i\Xi j}$ in the enclosing cell are used for the interpolation

$h_{i}(\mathbf{x},t)=\sum_{j=1}^{N_{j}}\hat{h}_{i\Xi j}(t)\psi_{\Xi
j}(\mathbf{x})\vspace{-.1cm}$ (8)

in cell $\Xi$ and with the basis functions $\psi_{\Xi j}$. A three-dimensional
reference cell with polynomial order $p=4$ is shown in Fig. 1.

Figure 1: Support points of a three-dimensional reference cell with Gauß-
Lobatto-Chebyshev points of order $p=4$. Interior and rear points are not
shown.

For each of the $N$ support points in the simulation, there are $Q$ particle
velocities, i.e. there are $N\cdot Q$ departure points to be identified.
Therefore, at the beginning of the simulation the path from each support point
to the corresponding $Q$ departure points is tracked through all adjacent
cells. Then a sparse matrix $\bm{\Psi}$ stores the interpolation coefficients
$\psi_{\Xi j}$ belonging to the departure point’s position; the algorithm is
presented in [43]. The actual streaming step is expressed as a matrix-vector
multiplication

$\bm{h}_{i}(t+\delta_{t})=\bm{\Psi}_{i}\bm{h}_{i}(t),$ (9)

whereas the collision step remains local.

###  II.2 Cubature-based velocity sets

The discretization of the velocity space is a key principle for any simulation
with the lattice Boltzmann method. If a quadrature is to be applied, it must
be suited to integrate the weight function $\mathrm{exp}{(-x^{2})}$, and it
has to be of ninth degree to enable compressible flow simulations [81]. A
prominent method to derive two- and three-dimensional velocity sets is the
Gauß-product rule applied to a one-dimensional quadrature. Application to the
one-dimensional degree-nine Gauß-Hermite quadrature delivers a two-dimensional
D2Q25 off-lattice velocity set with 25 abscissae, which we used for previous
work [44]. Due to its structure, this velocity set is infeasible for standard
on-lattice streaming but perfectly suited for the SLLBM.

For the simulations in this work we used velocity sets derived by cubature
rules [45, 46], exhibiting the same degree of precision but consisting of
fewer support points to lower computational cost. In two dimensions we
employed the degree-nine D2Q19 velocity set with 19 abscissae [82] that we
have presented in recent work [47]. In three dimensions, the Gauß-product rule
led to noncompetitive 125 abscissae for a three-dimensional degree-nine D3Q125
velocity set. Therefore, we successfully identified and implemented a degree-
nine D3Q45 velocity set with only 45 components, which fulfills all
requirements in terms of symmetry. The D3Q45 was derived by a cubature rule
after Konyaev [83] and has recently been listed in the supplemental material
of a cubature article by van Zandt [84]. The resulting discrete velocities are
shown in Fig. 2, whereas weights and abscissae of the D2Q19 and D3Q45 are
listed in [47].

Figure 2: Three-dimensional D3Q45 velocity set with 45 abscissae, derived by a
degree-nine cubature rule originally by Konyaev [83]

##  III Results

We test the proposed method through three test cases. The first two test
cases, a temporally underresolved Sod shock tube and a two-dimensional Riemann
problem, demonstrate the effect of extraordinarily large time step sizes on
the simulations. The third test case is the compressible three-dimensional
Taylor-Green vortex at different Reynolds and Mach numbers to present the
capability of SLLBM to simulate both turbulent and shocked three-dimensional
flows.

###  III.1 Temporally underresolved Sod shock tube

The Sod shock tube illustrates the large range of time step sizes accessible
to the SLLBM. The domain $x\in[0,1]$ was divided into two regions at $x=0.5$
with initially

$\displaystyle\rho_{0}=8,u_{0}=0,P_{0}=10,$
$\displaystyle\rho_{1}=1,u_{1}=0,P_{1}=1,$

and a viscosity of $\mu=7\cdot 10^{-4}$.

Figure 3: Density a) and velocity b) of Sod shock tube simulation using SLLBM
with 2000 cells at time $t=0.1$, which is reached after only 100 time steps
($\delta_{t}=0.001$), despite the fine spatial resolution of 6000 grid points.
For a smaller time step size and lower viscosity the reference is perfectly
matched. Figure 4: Exemplary departure point location for the Sod shock tube
configuration in this article with 2000 cells. The path of the longest
abscissa of the 2D velocity set D2Q19 is shown. Starting from the current
cell, the abscissa’s path linearly traverses six contiguous cells and locates
the departure point in the seventh cell. A 3rd-order polynomial interpolation
using the cell-local Gauß-Lobatto-Chebyshev support points is applied to
reconstruct the distribution function value.

The domain $x\in[0,1]$ was discretized using 2000 cells at third polynomial
order, i.e. 6000 grid points in x direction. Despite this fine spatial
resolution, which we only chose for demonstration purposes, the time step size
was set to $\delta_{t}=0.001$, such that the solution at $\mathrm{t}=0.1$ was
reached by performing only 100 time steps. Fig. 3 shows that the SLLBM
accurately traces the shock front despite the extremely large time step. For a
smaller time step size of $\delta_{t}=0.0005$ and viscosity of $\mu=10^{-5}$
the simulation results perfectly matched the inviscid reference solution. To
get an impression of the time step size of other solvers, we repeated the
simulations with the lower viscosity $\mu=7\cdot 10^{-4}$ using the same
collision process and velocity discretization, but this time applying a
spectral-element discontinuous Galerkin solver for the streaming step [74,
43]. This solver also features high-order solutions, but requires a dedicated
time integrator and the time step size is bounded by the CFL condition (1).
The simulation produced nearly identical results as the SLLBM with
$\delta_{t}=0.001$ (therefore not shown), but required $\delta_{t}=0.000086$
to be numerically stable, i.e., 1154 time steps with an explicit exponential
time integrator [85] and $\delta_{t}=0.00005$, i.e., 2200 time steps with the
more common fourth-order Runge-Kutta method.

As the distance of the departure points from the active nodes is proportional
to the time step size (Eq. 5), the departure points of the SLLBM were located
up to seven cells away. The trajectory is shown in Fig. 4 for an exemplary
departure point. It is obvious that the CFL restriction of explicit Eulerian
solvers prohibits the exchange of information crossing multiple cells. This
property is of special interest in the case of simulations with body-fitted
meshes, where the spatial extent of the smallest cells usually dictates the
time step size of the whole simulation. As opposed to Eulerian solvers, the
maximum stable time step size in the SLLBM is proportional only to the
viscosity and not dictated by the mesh size. On top of that, when doubling the
number of cells the number of SLLBM time steps can still be kept constant,
whereas it inevitably doubles for the explicit discontinuous Galerkin solver.
Finally, Fig. 5 confirms that the SLLBM is also capable to stably simulate the
shock tube with a lower resolution of 100 cells, polynomial order $p=4$ and
400 grid points, without additional numerical dissipation measures. Here, the
time step size was set to $\delta_{t}=0.005$ and $\delta_{t}=0.0005$ with
viscosities $\mu=0.002$ and $\mu=0.0002$, respectively.

Figure 5: a) Density b) velocity of Sod shock tube simulation using SLLBM with
100 cells at time $t=0.1$, which is reached after 20 time steps
($\delta_{t}=0.005$), at a more reasonable resolution of 400 grid points.
Simulations with smaller time step sizes and viscosity match the reference
better, despite slight overshoots visible near the shock fronts.

###  III.2 2D Riemann study of time step size effects

Case 12 of the two-dimensional Riemann problems was intensively studied by Lax
[86] as well as Kurganov and Tadmor [87]. In one of our last works [44] we
already showed that the SLLBM is capable to resolve the density contours of
this test case with good visual agreement to the references. This time,
similar to the shock tube test case in Section 3, we examined the effect of
the time step size onto this test case. The initial conditions of the four
quadrants are [76]

$\displaystyle(\rho,u_{x},u_{y},p)=\begin{cases}(0.5313,0,0,0.4),\;&x>0,\;y>0,\\\
(1,0.7276,0,1),\;&x\leq 0,\;y>0,\\\ (0.8,0,0,1),\;&x\leq 0,\;y\leq 0,\\\
(1,0,0.7276,1),\;&x>0,\;y\leq 0.\end{cases}$ (10)

The heat capacity ratio was $\gamma=1.4$, the number of cells was
$N_{\Xi}=128^{2}$ with polynomial order $p=4$, i.e., $N=512^{2}$ grid points.
The simulation end time was $t_{\mathrm{end}}=0.25$. Fig. 6 shows the density
contours in the interval $\rho\in[0.412,1.753]$, beginning with a small time
step size and low viscosity in simulation a). Each of the figures also lists
the time step size and viscosity. When increasing the time step size, the
viscosity also needs to be increased to ensure stable simulations as depicted
for simulation b). The effect of a maladjusted time step size is shown in
simulation c), where oscillations occur in the top right corner near the shock
fronts. Finally, simulation d) depicts the density contours for a
approximately 50 times larger time step size than for case a). In this case,
the shock fronts are widened, comparable to the observations in Section III.1.

Figure 6: Density contours of the 2D Riemann simulations using $512\times 512$
grid points. The time step size of simulation d) is approximately 50 times
larger than the time step size of simulation a), but the viscosity differs by
a factor of approximately 12. As a general rule, lower viscosities demand
smaller time step sizes. When choosing the time step size too large for a
given viscosity, the simulations tend to become unstable near the shock
fronts, as depicted in c). By contrast, simulation b) remains stable due to a
57 percent reduction of the time step size with equal prescribed viscosity.

###  III.3 Compressible Three-Dimensional Taylor-Green vortex

To show that the SLLBM captures the intricate interactions between turbulent
and compressible features including shocklets, the compressible three-
dimensional Taylor-Green vortex flow was simulated. On the triply-periodic
domain $S=[0,2\pi]^{3}$, the initial conditions with constant temperature
initial condition (CTIC) are

$\displaystyle u_{1}(x_{1},x_{2},x_{3},t=0)$ $\displaystyle=\sin x_{1}\cos
x_{2}\cos x_{3},$ (11) $\displaystyle u_{2}(x_{1},x_{2},x_{3},t=0)$
$\displaystyle=-\cos x_{1}\sin x_{2}\cos x_{3},$ (12) $\displaystyle
u_{3}(x_{1},x_{2},x_{3},t=0)$ $\displaystyle=0,$ (13)
$\displaystyle\rho(x_{1},x_{2},x_{3},t=0)$
$\displaystyle=1+\frac{\mathcal{C}}{16}[\cos(2x_{1})$ (14)
$\displaystyle+\cos(2x_{2})]\left[\cos(2x_{3}+2)\right],$ (15) $\displaystyle
T(x_{1},x_{2},x_{3},t=0)$ $\displaystyle=1,$ (16)

with velocities $\mathbf{u}$, density $\rho$, Mach number $\mathrm{Ma}$, and
temperature $T$. The numerator $\mathcal{C}$ differs between the cases in this
work with $\mathcal{C}=1$ for $\mathrm{Re}=100$ and $\mathrm{Re}=400$ as well
as $\mathcal{C}=\gamma\mathrm{Ma}^{2}$ for $\mathrm{Re}=1600$. The Reynolds
number is defined as $Re=1/\nu_{\infty}$, where the subscript $\infty$ denotes
the value at $T=1$. The Prandtl number is $\mathrm{Pr}=0.7$; the heat capacity
ratio is $\gamma=1.4$. The dynamic viscosity $\mu=\nu\rho$ obeys the
Sutherland law

$\mu=\frac{1.4042T^{1.5}}{T+0.40417}\mu_{\infty}.$ (17)

In comparison to forced or decaying isotropic turbulence, this test case
enables a deterministic initialization and thus an easier and more objective
comparison.

####  III.3.1 Reynolds number $\mathrm{Re}=400$

Peng and Yang thoroughly studied the compressible Taylor-Green vortex at
Reynolds number $\mathrm{Re}=400$ [88]. The original work used a compact
eighth-order finite difference scheme [24] to discretize the convective terms
in the Navier-Stokes equation in combination with a seventh-order weighted
essentially non-oscillatory (WENO) scheme. The present compressible SLLBM uses
fourth-order polynomials for the equilibrium (Eq. (3)) and for the
interpolation (Eq. (8)), _without utilizing additional filtering or
stabilizing_ techniques. Moreover, to satisfy the CFL condition, the original
work applied a time step size of $\delta_{t}=0.0005$, wheares the SLLBM was
capable to utilize time step sizes two orders of magnitude larger:
$\delta_{t}=0.017$ for Mach number $\mathrm{Ma}=0.5$, $\delta_{t}=0.033$ for
$\mathrm{Ma}=1.0$, $\delta_{t}=0.018$ for $\mathrm{Ma}=1.5$, and
$\delta_{t}=0.012$ for $\mathrm{Ma}=2.0$. The spatial resolution was
$N_{\mathrm{points}}=256^{3}$, whereas the reference operated with $512^{3}$
grid points.

Figure 7: Kinetic energy over time of the compressible 3D Taylor-Green vortex
at Reynolds number $\mathrm{Re}=400$. Simulations by the SLLBM with the D3Q45
velocity set shown in Fig. 2. Reference from [88].

Fig. 7 depicts the kinetic energy over time for all Mach numbers, which
accurately follows the reference solution despite the coarse temporal
discretization. Next, Figure 8 depicts the solenoidal dissipation defined as

$\varepsilon^{s}=\frac{1}{\rho_{\mathrm{ref}}\mathrm{Re}}\int_{S}\mu(\nabla\times\mathbf{u})^{2}\
d^{3}\mathbf{x}.$ (18)

All Mach numbers quite accurately matched the reference. The effect of
changing the time step to similar levels as in the above mentioned reference
can be seen for the smaller Mach numbers $\mathrm{Ma}=0.5$ and
$\mathrm{Ma}=1.0$. The small deviation from the reference seen here is reduced
by decreasing the time step size to similar levels, here $\delta_{t}=0.003$ in
both cases.

The dilatational dissipation

$\varepsilon^{c}=\frac{4}{3\rho_{\mathrm{ref}}\mathrm{Re}}\int_{S}\mu(\nabla\cdot\mathbf{u})^{2}\
d^{3}\mathbf{x}$ (19)

is a measure for pressure work in the simulation. Fig. 9 shows that the
dilatational effects are strong in the beginning at small simulation times,
surmounting the solenoidal dissipation. The comparison to the reference shows
a slight deviation from the reference and mild oscillations near the peak
values around $t=2.5$.

Figure 8: Solenoidal dissipation $\varepsilon^{s}$ of the compressible 3D
Taylor-Green vortex at Reynolds number $\mathrm{Re=400}$. Small time step
sizes at $\mathrm{Ma}=0.5$ and $\mathrm{Ma}=1.0$ with 5.5 times and 11.1 times
smaller times step sizes, respectively. The reduction induced a slightly
better agreement with the reference from [88]. Figure 9: Dilatational
dissipation $\varepsilon^{c}$ of the compressible 3D Taylor-Green vortex at
Reynolds number $\mathrm{Re=400}$. The time step sizes are identical to the
large $\delta_{t}$ simulations in Fig. 8.

####  III.3.2 Reynolds number $Re=1600$

The Taylor-Green vortex at Reynolds number $\mathrm{Re}~{}=~{}400$ shows
strong dilatational effects, but the transition to fully developed turbulence
requires higher Reynolds numbers. Therefore, another recent work by Lusher and
Sandham examined this test case at Reynolds number $\mathrm{Re}=1600$ up to
Mach number $\mathrm{Ma}=1.25$ [89]. In their study, the authors compared
high-order finite difference schemes equipped by WENO or targeted essentially
non-oscillatory (TENO) schemes of different orders. The present work re-
examines the Mach numbers $\mathrm{Ma}=1.0$ and $\mathrm{Ma}=1.25$ up to
$t=20$, with time step sizes $\delta_{t}=0.010$, and $\delta_{t}=0.009$,
respectively. These time step sizes were still 20 times larger than those in
the reference [89]. The resolution was $N_{\mathrm{points}}=256^{3}$, whereas
the reference used $512^{3}$ grid points. Fig. 10 demonstrates that, just as
in the case $\mathrm{Re}=400$, the kinetic energy over time is excellently
reproduced for both Mach numbers.

Figure 10: Kinetic energy over time of the compressible 3D Taylor-Green vortex
at Reynolds number $\mathrm{Re}=1600$ for Mach numbers a) $\mathrm{Ma}=1.0$
and b) $\mathrm{Ma}=1.25$. Reference from [89]. Figure 11: Solenoidal
dissipation $\varepsilon^{s}$ of the compressible 3D Taylor-Green vortex at
Reynolds number $\mathrm{Re=1600}$ for Mach numbers a) $\mathrm{Ma}=1.0$ and
b) $\mathrm{Ma}=1.25$. Note the differences between second-order polynomials
$p=2$ and fourth-order polynomials $p=4$.

Fig. 11 shows the solenoidal dissipation in comparison to the reference with
even better results for the higher Mach number $\mathrm{Ma}=1.25$. The good
agreement confirms the little numerical dissipation introduced by the SLLBM at
polynomial order $p=4$ in turbulent flows at transonic Mach numbers. By
contrast, the interpolation order $p=2$ worsens the solution despite the
identical resolution of $N_{\mathrm{points}}=256^{3}$.

Figure 12: Dilatational dissipation $\varepsilon^{c}$ of the compressible 3D
Taylor-Green vortex at Reynolds number $\mathrm{Re=1600}$.

Next, we investigated the dilatational dissipation. During the early phase of
the simulation the large vortex structures begin to entangle. This
entanglement leads to strong moving shock-like structures or turbulent
shocklets with strong negative dilatation [89] and local numerical
oscillations of the macroscopic variables due to an underresolution of these
shocklets. A slice of the density field illustrating the shocklets, indicated
by large jumps in the governing variables, is shown in Fig. 13. The size of
the jumps agreed well with that obtained via classical Rankine-Hugoniot jump
conditions. Additionally, during the early phase, the Mach numbers can be
higher than the initially prescribed Mach numbers. Moreover, to compute Eq.
(19), we made use of the gradients of the interpolation polynomials, which—in
contrast to the distribution functions—are not continuous over the cells. It
is likely that this approach to compute the dilatational dissipation faces
issues at strong shocklets. This explains the oscillations for
$\mathrm{Ma}=1.25$ in Fig. 12, which depicts the dilatational dissipation.
Despite these deviations from the reference in the beginning, the SLLBM was
able to reproduce the dilatational dissipation well for the rest of the
simulation, as shown in Fig. 12, without using additional stabilizing measures
like filtering or artificial diffusivity.

Figure 13: Density of the three-dimensional Taylor-Green vortex with
$\mathrm{Re}=1600$ and Mach number $\mathrm{Ma}=1.25$ at time $t~{}=~{}3.34$.
The x-z plane is shown at position $y=\pi$.

####  III.3.3 Reynolds number $\mathrm{Re}=100$ for higher Mach numbers

As a last test case, Fig. 14 shows that for $\mathrm{Re}=100$ we were able to
stably simulate the given test case for Mach numbers $\mathrm{Ma}=2.5$ and
$\mathrm{Ma}=3.0$. The Knudsen number for this configuration was of order
$\mathrm{Kn}=\mathcal{O}(\mathrm{Ma}/\mathrm{Re})=\mathcal{O}(10^{-2})$. Due
to this rather large Knudsen number, the SLLBM was once more able to use large
times steps $\delta_{t}=0.015$ and $\delta_{t}=0.03$, respectively. Although
only shown at low Reynolds numbers, these stable simulations indicate the
principal capability of the SLLBM to perform stable simulations at high Mach
numbers.

Figure 14: Kinetic energy over time of the compressible 3D Taylor-Green vortex
at Reynolds number $\mathrm{Re}=100$ for Mach numbers $\mathrm{Ma}=2.5$ and
$\mathrm{Ma}=3.0$.

##  IV Discussion

As demonstrated by the numerical experiments, the SLLBM is able to simulate
complex two- and three-dimensional compressible flows with shocks and
shocklets. There are two main arguments in favor of the method. First, the
method works with large and adjustable time step sizes. Second, cubature rules
provide lean velocity sets resulting in an efficient scheme. The D3Q45 used in
this work is, to the authors knowledge, the smallest known degree-nine
velocity set and key to reduce errors caused by interpolation.

The SLLBM is not subject to a CFL condition, as already shown in past works
[42] and as confirmed by the shock tube simulations in this work. In contrast
to the usual spatial filtering, the independence of the CFL condition opens
the field of temporal large eddy simulations (TLES) with fine spatial, but
coarse temporal resolutions [90, 91]. Since no dedicated filter operation has
been used, we classify the scheme as some sort of temporal implicit large eddy
simulation (TILES). Despite the temporally coarse resolutions, most relevant
flow features during the simulations were preserved, or they were recovered by
scaling the time step size down: a key feature of the SLLBM that on-lattice
Boltzmann schemes often miss.

Although there is no CFL condition, the simulations are not automatically
stable, because the collision step may still be prone to instabilities at
small relaxation times as pointed out by various works, e.g. by Siebert et al.
[92]. That said, the stability regions of the BGK collision operator combined
with the D3Q45 and fourth-order equilibria (at high Mach numbers and
adjustable time step sizes) are yet unknown. Nevertheless, we observed that
the SLLBM produces good results without using additional filtering or
stabilizing techniques. To accomplish simulations with even higher Reynolds
and Mach numbers with large time step sizes, the application of such
techniques to the SLLBM might be beneficial, though.

###  IV.1 Comparison to other LBM solvers

Let us now compare the SLLBM for three-dimensional compressible flows to other
LBM solvers, with distinction between on-lattice solvers and off-lattice
solvers. Since we discussed many aspects in recent works [44, 47], we mainly
restrict the discussion to compressible LBM solvers with applications to
three-dimensional turbulent flows.

To perform a time step, off-lattice Boltzmann methods typically require a
special treatment of the distribution function values such as the
discretization by finite difference schemes [93, 94, 95] or finite volume
schemes [76] and the application of a dedicated time integrator, e.g., a
Runge-Kutta scheme. As an example for these Eulerian time integration schemes,
Chen et al. presented compressible decaying three-dimensional isotropic
turbulence simulations obtained by a modified discrete unified gas kinetic
scheme (DUGKS), which is essentially a finite volume LBM with second-order
spatial and temporal accuracy [96]. Also, the authors made use of a high-order
Gaussian quadrature to discretize the velocity space, but they relied on a
decisively larger D3Q77 velocity set with identical quadrature degree as the
D3Q45. Like all explicit Eulerian time integration schemes, DUGKS is subject
to the CFL condition (1).

A second category of off-lattice schemes are interpolation-based or semi-
Lagrangian implementations. Pavlo et al. pioneered in using second-order
interpolations to incorporate high-order velocity sets for simulations of
thermal flows [97, 98]. Renowned semi-Lagrangian implementations are the
Particles on Demand method [61], which uses dynamically shifted velocity sets
to reach high Mach numbers or the SLLBM for compressible flows on unstructured
meshes by Saadat et al. [65]. Unlike the present method, the authors of the
latter approach used a D2Q9 velocity set and computed correction terms by
exploiting the gradients of the distribution functions, which practically come
at low costs when using interpolation polynomials. To the best of our
knowledge, these methods have not been applied to complex three-dimensional
flows, except for a spherical Riemann problem presented by Zakirov et al.
using a D3Q125 velocity set [99]. Additionally, one decisive feature of the
present method is the organization of support points by cells, enabling three-
dimensional high-order solutions.

On-lattice Boltzmann solvers have the Lagrangian integration along
characteristics in common with interpolation-based schemes, albeit constrained
by a strong coupling of space and time discretization. This explains, the
rather large time step sizes of on-lattice Boltzmann methods with increasing
Mach number [100, 47]. On-lattice Boltzmann methods generally exhibit second-
order accuracy in space and time, but they are not as fiercely concerned by
numerical diffusion as low-order off-lattice Boltzmann methods. In general,
the computational costs of on-lattice Boltzmann methods’ streaming step are
low. On the downside, they often suffer from large velocity sets. For example,
Frapolli et al. [55, 57, 101] were able to simulate various three-dimensional
compressible flows, including isotropic decaying turbulence simulations or the
flow around an Onera M6 wing, by the entropic lattice Boltzmann method. In
many ways, their works were groundbreaking for the development of compressible
LBMs. However, the authors used a velocity set with 343 discrete velocities
and the weights were temperature-dependent, which limits the temperature range
of the method. In addition, in many on-lattice Boltzmann frameworks the time
step size is not adjustable due to the linking of spatial and temporal
discretization [100]. An exception to this are hybrid lattice Boltzmann
methods, which solve the temperature field separately [102, 58, 62, 103]. In
this case, by adapting the speed of sound, the time step size can be varied.
By using numerical equilibria instead of the polynomial equilibrium in this
work, Latt et al. were able to simulate a three-dimensional flow around a
sphere [66] by using only 39 discrete velocities, as similarly proposed by
Frapolli [100]. Simulations of complex turbulent flows by this method are not
available, yet. Recently, Saadat et al. [104] proposed an on-lattice Boltzmann
model with a regular D3Q27 velocity set to perform three-dimensional decaying
compressible isotropic turbulence simulations. Despite the low degree of the
velocity set, this model proved capable to simulate shocked flows up to Mach
numbers of $\mathrm{Ma}~{}1.5$. This is achieved by correcting the error-prone
high-order moments by expressions obtained from applying finite differences to
the macroscopic variables.

In summary, the number of compressible LBM solvers for three-dimensional flows
is still limited for both on-lattice and off-lattice Boltzmann schemes, which
indicates the value of the presented SLLBM framework to set a pattern for
future research.

###  IV.2 Numerical efficiency

The numerical efficiency of the SLLBM depends on the implementation of the
collision and the streaming step. The former requires at each support point
the calculation of the discrete equilibrium function values given by Eq. (6).
Once the density, velocities, and temperature are gained from the distribution
function values, the determination of the equilibrium, however, is well
parallelizable. Note, that the tensors in Eq. (6) are symmetric, so that many
entries need to be computed only once per time step and support point, e.g.,
$a_{xxxy,\mathrm{eq}}$ is equal to $a_{yxxx,\mathrm{eq}}$ and all other index
permutations. In addition, the entries of the Hermite tensors
$\mathcal{H}_{i}$ are constant for a given velocity set and only need to be
determined once for the whole simulation.

Another decisive factor of the simulation performance is the streaming step,
which requires a “triquartic” interpolation, i.e., interpolation polynomials
of order $p=4$ in three dimensions. The reference cell is shown in Fig. 1. To
interpolate one distribution function value at a given departure point, all
$(p+1)^{D}=125$ support points of the cell are processed with the
interpolation coefficients that are stored in the $Q$ matrices
$\bm{\Psi}_{i}$. These matrices are the cornerstone for the whole simulation.
When accounting to equally shaped cells, like in all simulations in the
present work, the size of the matrices can be significantly reduced, since the
matrices are identical over nearly all cells (boundary cells excepted).
However, for irregular or distorted three-dimensional meshes, the size of the
matrices grows quickly, rendering the streaming step memory-bound. Therefore,
matrix-free implementations [105] appear to be an attractive extension for
SLLBM implementations. In matrix-free versions the interpolation coefficients
will be computed afresh in each time step, which is potentially faster for
today’s high-performance computing clusters.

At first glance, the second distribution function appears to be a severe
limitation of the method. On closer inspection, however, it turns out that the
overhead is kept in reasonable limits since the computation-intensive
equilibrium distribution $g_{i}^{\mathrm{eq}}$ linearly depends on
$f_{i}^{\mathrm{eq}}$ and the interpolation coefficients can be used for both
$f_{i}$ and $g_{i}$. Still, all aforementioned issues pay off in light of the
large time step sizes of semi-Lagrangian implementations. A major way to
further reduce the computational cost is the reduction of discrete velocities
in the Gaussian quadrature. The research for even compacter cubature rules is
still ongoing [47], so that future degree-nine velocity sets will possibly be
even more efficient.

One last remark regarding the equilibrium distribution function determined via
Eq. (6). In the past the Hermite expansion-based equilibrium was deemed
causing instabilities due to negative distribution function values at high
Mach numbers [100]. Fig. 15 depicts the equilibrium function values
$f^{\mathrm{eq}}_{i}$ over the Mach number $\mathrm{Ma}_{x}$ in x-direction.
The figure manifests that most of the discrete velocities significantly
diverge for Mach numbers larger than $\mathrm{Ma}>2.0$ towards either $\infty$
or $-\infty$. However, our simulations showed that simulations even up to
$\mathrm{Ma}=3.0$ remained stable. This observation indicates that negative
distribution function values are not per se a source of instabilities, since
the values’ sole role is to encode the moments of different orders. For the
first moments this encoding can also be seen in Fig. 15: despite the
increasing first-order moment of the momentum, the “zeroth–order” moment of
the density remains $\rho=1.0$. To ensure this, negative values are a
consequence at high Mach numbers due to the discretization of the velocity
space via Gaussian quadratures.

Figure 15: All 45 equilibrium distribution functions $f^{\mathrm{eq}}_{i}$ of
the D3Q45 velocity set over Mach numbers ranging from 0 to 4 in x-direction.
Density $\rho$ and relative temperature $\theta$ were constantly set to unity;
Mach numbers in both other directions were set to zero.

##  V Conclusion

The SLLBM for three-dimensional compressible flows is a viable alternative to
other solvers. The SLLBM allows very large time steps sizes not restricted by
the customary CFL condition. Although the presented SLLBM is a fourth-order
spatial method and accurately captures shocks as well as turbulence, no
stabilization or filtering were required for the presented test cases. Due to
these unique features, the cubature-based fully explicit SLLBM enables
researchers to perform compressible turbulence simulations, in which the
admissible time step sizes are decoupled from the spatial discretization,
opening a new field of affordable simulations for compressible turbulent
flows.

###### Acknowledgements.

This work was supported by the German Ministry of Education and Research and
the Ministry for Culture and Science North Rhine-Westfalia (research grant
13FH156IN6). D.W. is supported by German Research Foundation (DFG) project FO
674/17-1.
D.W. and A.K. contributed equally to this work.

##  Appendix A Equilibrium moments

The equilibrium moments up to fourth order read

$\displaystyle a^{(0)}_{eq}$ $\displaystyle=$ $\displaystyle\rho$ (20a)
$\displaystyle a^{(1)}_{\alpha,eq}$ $\displaystyle=$ $\displaystyle\rho
u_{\alpha}$ (20b) $\displaystyle a^{(2)}_{\alpha\beta,eq}$ $\displaystyle=$
$\displaystyle\Pi_{\alpha\beta}^{\mathrm{eq}}=\rho(u_{\alpha}u_{\beta}+T_{0}(\theta-1)\delta_{\alpha\beta})$
(20c) $\displaystyle a^{(3)}_{\alpha\beta\gamma,eq}$ $\displaystyle=$
$\displaystyle\mathcal{Q}_{\alpha\beta\gamma}^{\mathrm{eq}}=\rho\left[u_{\alpha}u_{\beta}u_{\gamma}+T_{0}(\theta-1)(\delta_{\alpha\beta}u_{\gamma}\right.$
(20d) $\displaystyle\left.\hskip
34.14322pt+\delta_{\alpha\gamma}u_{\beta}+\delta_{\beta\gamma}u_{\alpha})\right]$
$\displaystyle a^{(4)}_{\alpha\beta\gamma\delta,eq}$ $\displaystyle=$
$\displaystyle\mathcal{R}_{\alpha\beta\gamma\delta}^{\mathrm{eq}}$ (20e)
$\displaystyle=$
$\displaystyle\rho[u_{\alpha}u_{\beta}u_{\gamma}u_{\delta}+T_{0}(\theta-1)((\delta_{\alpha\beta}\delta_{\gamma\delta}$
$\displaystyle+\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma})T_{0}(\theta-1)$
$\displaystyle+\delta_{\alpha\beta}u_{\gamma}u_{\delta}+\delta_{\alpha\gamma}u_{\beta}u_{\delta}+\delta_{\alpha\delta}u_{\beta}u_{\gamma}+$
$\displaystyle\delta_{\beta\gamma}u_{\alpha}u_{\delta}+\delta_{\beta\delta}u_{\alpha}u_{\gamma}+\delta_{\gamma\delta}u_{\alpha}u_{\beta})],\vspace{.3cm}$

where $\theta=T/T_{0}$ is the relative temperature.

##  Appendix B Hermite tensors

The scaled Hermite tensors with $\hat{\xi_{i}}=\xi_{i}/c_{s}$ up to fourth
order read

$\displaystyle{\mathcal{H}_{i}}^{(0)}$ $\displaystyle=1$
$\displaystyle{\mathcal{H}_{i\alpha}}^{(1)}$
$\displaystyle=\frac{{{\hat{\xi}_{i\alpha}}}}{c_{s}}$
$\displaystyle{\mathcal{H}_{i\alpha\beta}}^{(2)}$
$\displaystyle=\frac{{{{\hat{\xi}_{i\alpha}}{\hat{\xi}_{i\beta}}}}-\delta_{\alpha\beta}}{c_{s}^{2}}$
$\displaystyle{\mathcal{H}_{i\alpha\beta\gamma}}^{(3)}$
$\displaystyle=\frac{{\hat{\xi}_{i\alpha}}{\hat{\xi}_{i\beta}}{\hat{\xi}_{i\gamma}}-({\hat{\xi}_{i\alpha}\delta_{\beta\gamma}}+{\hat{\xi}_{i\beta}\delta_{\alpha\gamma}}+{\hat{\xi}_{i\gamma}\delta_{\alpha\beta}})}{c_{s}^{3}}$
$\displaystyle{\mathcal{H}_{i\alpha\beta\gamma\delta}}^{(4)}$
$\displaystyle=\frac{{{\hat{\xi}_{i\alpha}}{\hat{\xi}_{i\beta}}{\hat{\xi}_{i\gamma}}{\hat{\xi}_{i\delta}}}-\mathcal{T}_{i}+(\delta_{\alpha\beta}\delta_{\gamma\delta}+\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma})}{c_{s}^{4}},$

with

$\begin{split}\mathcal{T}_{i}=\hat{\xi}_{i\alpha}\hat{\xi}_{i\beta}\delta_{\gamma\delta}+\hat{\xi}_{i\alpha}\hat{\xi}_{i\gamma}\delta_{\beta\delta}+\hat{\xi}_{i\alpha}\hat{\xi}_{i\delta}\delta_{\beta\gamma}+\\\
\hat{\xi}_{i\beta}\hat{\xi}_{i\gamma}\delta_{\alpha\delta}+\hat{\xi}_{i\beta}\hat{\xi}_{i\delta}\delta_{\alpha\gamma}+\hat{\xi}_{i\gamma}\hat{\xi}_{i\delta}\delta_{\alpha\beta}.\end{split}$

##  Appendix C Chapman-Enskog analysis

This section shows the approximation of the compressible Navier-Stokes
equations, when applying a Chapman-Enskog analysis to the SLLBM model. To that
end, a second-order multiscale expansion is applied to the temporal
derivatives

$\displaystyle\partial_{\alpha}$ $\displaystyle=\epsilon\partial_{\alpha},$
(21) $\displaystyle\partial_{t}$
$\displaystyle=\epsilon\delta_{t}^{(1)}+\epsilon^{2}\delta_{t}^{(2)},$ (22)

where $\epsilon$ is a smallness parameter usually identified as the Knudsen
number [79] and the discrete distribution functions $h_{i}$

$h_{i}=h_{i}^{(0)}+\epsilon h_{i}^{(1)}+\epsilon^{2}h_{i}^{(2)}.$ (23)

These expansion terms are applied to a second-order Taylor expansion of Eq.
(5)

$(\epsilon
D_{i}+\frac{\epsilon^{2}}{2}D_{i}D_{i})h_{i}=-\frac{1}{\tau}\left(h_{i}-h_{i}^{\mathrm{eq}}\right),$
(24)

with the material derivative [106]

$\epsilon D_{i}=\epsilon
D_{i}^{(1)}+\epsilon^{2}D_{i}^{(2)}=\epsilon\partial_{t}^{(1)}+\epsilon\,\xi_{i\alpha}\partial_{\alpha}+\epsilon^{2}\partial_{t}^{(2)},$
(25)

Next, Eq. (23) is applied to (24) and the terms of same order are collected.
The zeroth-order terms of order $\epsilon^{0}$ are

$0=-\frac{1}{\tau}\left(h^{(0)}-h^{\mathrm{eq}}\right),$ (26)

with the trivial relation

$h_{i}^{(0)}=h_{i}^{\mathrm{eq}}.$ (27)

When applying the moment relations of Eq. (7) to the last Equation, we need to
differ between the moments of $f_{i}$ and $g_{i}$, i.e,

$\sum_{i}f_{i}^{(0)}\\{1,\xi_{i\alpha}\\}=\sum_{i}f_{i}^{\mathrm{eq}}\\{1,\xi_{i\alpha}\\}$
(28)

and

$\sum_{i}\left(f_{i}^{(0)}|\bm{\xi}_{i}|^{2}+g_{i}^{(0)}\right)=\sum_{i}\left(f_{i}^{\mathrm{eq}}|\bm{\xi}_{i}|^{2}+g_{i}^{\mathrm{eq}}\right),$
(29)

showing mass, momentum and energy conservation during the collision step,
which also implies that

$\displaystyle\sum_{i}f_{i}^{(1,2,\ldots)}\\{1,\xi_{i\alpha}\\}=0,$ (30)
$\displaystyle\sum_{i}\left(f_{i}^{(1,2,\ldots)}|\bm{\xi}_{i}|^{2}+g_{i}^{(1,2,\ldots)}\right)=0$
(31)

As a second step, the order $\epsilon^{1}$ terms are collected

$D_{i}^{(1)}h_{i}^{(0)}=-\frac{1}{\tau}h_{i}^{(1)}.$ (32)

This time, the relations of Eq. (7) yield for the zeroth moment

$D_{i}^{(1)}\rho=\partial_{t}^{(1)}\rho+\partial_{\alpha}(\rho u_{\alpha})=0,$
(33)

for the first moments with the help of Eq. (20c) and $P=\rho(T-T_{0})$

$\displaystyle\partial_{t}^{(1)}(\rho
u_{\alpha})+\partial_{\beta}\Pi^{(0)}_{\alpha\beta}=0$
$\displaystyle\partial_{t}^{(1)}(\rho u_{\alpha})+\partial_{\beta}(\rho
u_{\alpha}u_{\beta})+\partial_{\alpha}(P)\delta_{\alpha\beta}=0,$ (34)

and for the total energy

$\partial^{(1)}_{t}(\rho E)+\partial_{\alpha}(u_{\alpha}\rho
E)=-P\partial_{\alpha}u_{\alpha}$ (35)

or

$\partial^{(1)}_{t}T=-u_{\alpha}\partial_{\alpha}T-\frac{P}{\rho
C_{v}}\partial u_{\alpha}$ (36)

for the temperature.

These are the Euler equations; to derive the Navier-Stokes equations, the
terms of order $\epsilon^{2}$ also need to be gathered. This leads to

$D_{i}^{(2)}h_{i}^{(0)}+D_{i}^{(1)}h_{i}^{(1)}+\frac{1}{2}D_{i}^{(1)}D_{i}^{(1)}h_{i}^{(0)}=-\frac{1}{\tau}h_{i}^{(2)}$
(37)

Eq. (32) simplifies the last Equation’s derivatives

$D_{i}^{(2)}h_{i}^{(0)}+\left(1-\frac{1}{2\tau}\right)D_{i}^{(1)}h_{i}^{(1)}=-\frac{1}{\tau}h_{i}^{(2)}.$
(38)

Due to vanishing terms, the zeroth moment of Eq. (38) is

$\partial_{t}^{(2)}\rho=0.$ (39)

Then, the slightly more complicated part begins with the first moments of (38)

$\partial_{t}^{(2)}(\rho
u_{\alpha})=-\left(1-\frac{1}{2\tau}\right)\partial_{\beta}\Pi_{\alpha\beta}^{(1)}.$
(40)

To express the nonequilibrium moment $\Pi_{\alpha\beta}^{(1)}$ by equilibrium
moments, we derive the second order moment of (32), which is

$\partial_{t}^{(1)}\Pi_{\alpha\beta}^{(0)}+\partial_{\gamma}\mathcal{Q}_{\alpha\beta\gamma}^{(0)}=-\frac{1}{\tau}\Pi_{\alpha\beta}^{(1)}$
(41)

and then complements Eq. (40)

$\partial_{t}^{(2)}(\rho
u_{\alpha})=\left(\tau-\frac{1}{2}\right)\partial_{\beta}\left(\partial^{(1)}_{t}\Pi_{\alpha\beta}^{(0)}+\partial_{\gamma}\mathcal{Q}_{\alpha\beta\gamma}^{(0)}\right).$
(42)

The equilibrium moments can be explicitly specified by Equations (20c) and
(20d). By applying the product rule and by using the Euler equations for mass
(33) and momentum (34) and temperature (36) to replace the time derivatives,
one obtains

$\displaystyle\partial^{(1)}_{t}\Pi_{\alpha\beta}^{(0)}=$
$\displaystyle-\partial_{\gamma}(\rho
u_{\alpha}u_{\beta}u_{\gamma})-\frac{P}{C_{v}}\partial_{\gamma}u_{\gamma}\delta_{\alpha\beta}$
$\displaystyle-\partial_{\gamma}(Pu_{\gamma})\delta_{\alpha\beta}-u_{\alpha}\partial_{\beta}P-u_{\beta}\partial_{\alpha}P$
(43)

and

$\displaystyle\partial_{\gamma}\mathcal{Q}_{\alpha\beta\gamma}^{(0)}=$
$\displaystyle\partial_{\gamma}(\rho
u_{\alpha}u_{\beta}u_{\gamma})+P\partial_{\beta}u_{\alpha}+P\partial_{\alpha}u_{\beta}$
$\displaystyle+\partial_{\gamma}(Pu_{\gamma})\delta_{\alpha\beta}+u_{\alpha}\partial_{\beta}P+u_{\beta}\partial_{\alpha}P.$
(44)

This turns Eq. (42) into

$\partial_{t}^{(2)}(\rho u_{\alpha})=\\\
\left(\tau-\frac{1}{2}\right)\partial_{\beta}\left[P(\partial_{\alpha}u_{\beta}+\partial_{\beta}u_{\alpha})-\frac{P}{C_{v}}\partial_{\gamma}u_{\gamma}\delta_{\alpha\beta}\right].$
(45)

As a last step, the total energy is determined

$\partial_{t}^{(2)}(\rho
E)=-\left(1-\frac{1}{2\tau}\right)\partial_{\beta}q_{\beta}^{(1)},$ (46)

where

$q_{\beta}^{(1)}=\sum_{i}\frac{1}{2}|\bm{\xi}_{i}|^{2}\xi_{i\beta}f_{i}^{(1)}$
(47)

is a contracted variant of $\mathcal{Q}_{\alpha\beta\gamma}^{(1)}$ detailed in
Eq. (20d). Similar to Eq. (40), this vector can be expressed by equilibrium
counterparts

$q_{\beta}^{(1)}=-\tau(\partial_{t}^{(0)}q_{\beta}^{(0)}+\partial_{\gamma}r_{\beta\gamma}^{(0)}).$
(48)

with

$r_{\beta\gamma}^{(0)}=\sum_{i}\frac{1}{2}|\bm{\xi}_{i}|^{2}\xi_{i\beta}\xi_{i\gamma}f_{i}^{(0)}.$
(49)

Resembling $q_{\beta}$, the tensor $r_{\beta\gamma}$ is the contracted variant
of $\mathcal{R}_{\alpha\beta\gamma\delta}$ detailed in Eq. (20e). Again, by a
number of replacements, one obtains [107]

$q_{\beta}^{(1)}=-\tau P\bigl{[}(1+C_{v})\partial_{\beta}T\\\
+\partial_{\beta}\bigl{(}u_{\gamma}[\partial_{\beta}u_{\gamma}+\partial_{\gamma}u_{\beta}]-u_{\gamma}\frac{\partial_{\delta}u_{\delta}}{C_{v}}\delta_{\beta\gamma}\bigr{)}\bigr{]}.$
(50)

This term now complements Eq. (46). Finally, by summing up all contributions
of orders $\epsilon^{0}$, $\epsilon^{1}$, and $\epsilon^{2}$, the compressible
Navier-Stokes equations are derived

$\displaystyle\partial_{t}\rho+\partial_{\alpha}(\rho u_{\alpha})=0,$ (51)
$\displaystyle\partial_{t}(\rho u_{\alpha})+\partial_{\beta}(\rho
u_{\alpha}u_{\beta}+P\delta_{\alpha\beta})=\partial_{\beta}(\sigma_{\alpha\beta}),$
(52) $\displaystyle\partial_{t}(\rho E)+\partial(\rho
Eu_{\beta})=\partial_{\beta}(\lambda\partial_{\beta}T)+\partial_{\beta}(u_{\gamma}\sigma_{\beta\gamma}),$
(53)

with dynamic viscosity $\mu=\tau P$, thermal conductivity $\lambda=\tau
P(C_{v}+1)=\tau PC_{p}$, since the heat capacity at constant pressure is
defined as $C_{p}=C_{v}+1$. The stress tensor denotes

$\sigma_{\alpha\beta}=\mu\bigl{(}(\partial_{\alpha}u_{\beta}+\partial_{\beta}u_{\alpha})-P\delta_{\alpha\beta}\bigr{)}-\frac{\mu}{C_{v}}\partial_{\gamma}u_{\gamma}\delta_{\alpha\beta}.$
(54)

Since the dynamic viscosity depends on the local pressure $P$, the relaxation
parameter has to be pressure-dependent, i.e.,
$\tau=\mu/(Pc_{s}\delta_{t})+0.5$. From Eq. (54) the bulk viscosity, with
respect to the notation in Eq. (52)), can be identified as
$\mu_{b}~{}=~{}\mu/C_{v}$.

Note that the presented derivation links the thermal conductivity to the
dynamic viscosity. In our approach, we used the quasi-equilibrium approach to
adjust the Prandtl number $\mathrm{Pr}=(\mu C_{p})/\lambda$. For the
corresponding extension of the Chapman-Enskog analysis, we refer to [56].

##  Appendix D Quasi-equilibrium approach for variable Prandtl number

To obtain a variable Prandtl number $\mathrm{Pr}$, the following equation
replaces Eq. (5)

$h_{i}(\mathbf{x}\\!+\\!\delta_{t}\bm{\xi}_{i},t\\!+\\!\delta_{t})=h_{i}(x,t)\\\
-\frac{1}{\tau}\left[h_{i}(x,t)-h_{i}^{\mathrm{eq}}(x\xi_{i},t)\right]+\left(\frac{1}{\tau}-\frac{1}{\tau_{\mathrm{Pr}}}\right)h_{i}^{*}(x,t),$
(55)

with $\tau_{\mathrm{Pr}}=(\tau-0.5)/\mathrm{Pr}+0.5$. The quasi-nonequilibrium
$h_{i}^{*}$ is thereby obtained by first computing the centered heat flux
tensor

$\bar{Q}_{\alpha\beta\gamma}=\sum_{i=1}^{n}f_{i}(\xi_{i\alpha}-u_{\alpha})(\xi_{i\beta}-u_{\beta})(\xi_{i\gamma}-u_{\gamma}).$
(56)

Then, the nonequilibrium part
$\bar{Q}_{\alpha\beta\gamma}^{\mathrm{neq}}=\bar{Q}_{\alpha\beta\gamma}-\bar{Q}_{\alpha\beta\gamma}^{\mathrm{eq}}$
is applied to the Hermite tensor $\mathcal{H}_{i\alpha\beta\gamma}^{(3)}$

$h^{*}_{i}=w_{i}\bar{Q}_{\alpha\beta\gamma}^{\mathrm{neq}}:\mathcal{H}_{i\alpha\beta\gamma}^{(3)},$
(57)

via full contraction of indices.

## References

  * Lele [1994] S. K. Lele, Annual Review of Fluid Mechanics 26, 211 (1994).
  * Chassaing  _et al._ [2002] P. Chassaing, R. Antonia, F. Anselmet, L. Joly, and S. Sarkar,  _Variable Density Fluid Turbulence_ (Kluwer Academic Publishers, 2002).
  * Veynante [2002] L. Veynante, D. Vervisch, Prog. Energy Combust. Sci. 30, 193 (2002).
  * Smits and Dussauge [2006] A. J. Smits and J. P. Dussauge,  _Turbulent Shear Layers in Supersonic Flow: Second Edition_ (Springer, 2006).
  * Garnier  _et al._ [2009] E. Garnier, N. Adams, and P. Sagaut,  _LES for compressible flows_ (Springer, 2009).
  * Pirozzoli [2011] S. Pirozzoli, Annu. Rev. Fluid Mech. 43, 163 (2011).
  * Gatski and Bonnet [2013] T. B. Gatski and J.-P. Bonnet,  _Compressibility, Turbulence and High Speed Flow_ (Academic Press, 2013).
  * Sharman and Lane [2016] R. Sharman and T. Lane,  _Aviation turbulence: Processes, detection, prediction_ (Springer, 2016) pp. 1–523.
  * Low and Klessen [2004] M. M. M. Low and R. S. Klessen, Rev. Mod. Phys. 76, 125 (2004), 0301093 [astro-ph] .
  * Fernholz and Finley [1976] H. H. Fernholz and P. J. Finley,  _A critical compilation of compressible turbulent boundary layer data_ (AGARDograph 223, 1976).
  * Huang  _et al._ [1995] P. Huang, G. Coleman, and P. Bradshaw, J. Fluid Mech. 305, 185 (1995).
  * Foysi  _et al._ [2004] H. Foysi, S. Sarkar, and R. Friedrich, J. Fluid Mech. 509, 207 (2004).
  * Ghosh  _et al._ [2010] S. Ghosh, H. Foysi, and R. Friedrich, J. Fluid Mech. 648, 155 (2010).
  * Bradshaw [1977] P. Bradshaw, Ann. Rev. Fluid Mech. 21, 926 (1977).
  * Bogdanoff [1983] D. W. Bogdanoff, AIAA Journal 21, 926 (1983).
  * Vreman  _et al._ [1996] A. W. Vreman, N. D. Sandham, and K. H. Luo, J. Fluid Mech. 320, 235 (1996).
  * Foysi and Sarkar [2010] H. Foysi and S. Sarkar, Theoretical and Computational Fluid Dynamics 24, 565 (2010).
  * Matsuno and Lele [2020] K. V. Matsuno and S. K. Lele, AIAA 2020-0573, AIAA SciTech 2020 Forum (2020).
  * Colonius and Lele [2004] T. Colonius and S. Lele, Progr. Aerosp. Sc. 40, 345 (2004).
  * Bodony and Lele [2006] D. Bodony and S. Lele, AIAA, AIAA-2006-0468 (2006).
  * Foysi  _et al._ [2010] H. Foysi, M. Mellado, and S. Sarkar, International Journal of Heat and Fluid Flow 31, 307 (2010).
  * Lele and Larsson [2009] S. K. Lele and J. Larsson, J. Phys. Conf. Ser. 180, 012032 (2009).
  * Lee  _et al._ [1991] S. Lee, S. K. Lele, and P. Moin, Phys. Fluids A Fluid Dyn. 3, 657 (1991).
  * Lele [1992] S. K. Lele, J. Comput. Phys. 103, 16 (1992).
  * Tam and Webb [1993] C. K. Tam and J. C. Webb, J. Comput. Phys. 107, 262 (1993).
  * Adams and Shariff [1996] N. A. Adams and K. Shariff, J. Comput. Phys. 127(1) (1996).
  * Ekaterinaris [1999] J. A. Ekaterinaris, J. Comput. Phys. 156, 272 (1999).
  * Sengupta  _et al._ [2003] T. K. Sengupta, G. Ganeriwal, and S. De, J. Comput. Phys. 192, 677 (2003).
  * Bogey and Bailly [2004] C. Bogey and C. Bailly, J. Comput. Phys. 194, 194 (2004).
  * Wang  _et al._ [2010] J. Wang, L. P. Wang, Z. Xiao, Y. Shi, and S. Chen, J. Comput. Phys. 229, 5257 (2010).
  * Hu  _et al._ [1996] F. Hu, M. Hussaini, and J. Manthey, J. Comput. Phys. 124, 177 (1996).
  * Berland  _et al._ [2006] J. Berland, C. Bogey, and C. Bailly,  _Low-dissipation and low-dispersion fourth-order Runge-Kutta algorithm._ , Vol. 35 (Comput. Fluids, 2006) pp. 1459–1463.
  * Kwatra  _et al._ [2009] N. Kwatra, J. Su, J. T. Grétarsson, and R. Fedkiw, J. Comput. Phys. 10.1016/j.jcp.2009.02.027 (2009).
  * Müller [1990] B. Müller, Math. Models Methods Appl. Sci. 12, 139 (1990).
  * Xiu and Karniadakis [2001] D. Xiu and G. E. Karniadakis, J. Comput. Phys. 172, 658 (2001).
  * Celledoni  _et al._ [2016] E. Celledoni, B. K. Kometa, and O. Verdier, J. of Sci. Computing 66, 91 (2016), 1207.5147 .
  * Sonnendrücker  _et al._ [1999] E. Sonnendrücker, J. Roche, P. Bertrand, and A. Ghizzo, J. Comput. Phys. 149, 201 (1999).
  * Crouseilles  _et al._ [2010] N. Crouseilles, M. Mehrenberger, and E. Sonnendrücker, J. Comput. Phys. 229, 1927 (2010).
  * Qiu and Christlieb [2010] J. M. Qiu and A. Christlieb, J. Comput. Phys. 229, 1130 (2010).
  * Russo  _et al._ [2012] G. Russo, P. Santagati, and S. B. Yun, SIAM J. Num. Analysis 50, 1111 (2012).
  * Groppi  _et al._ [2014] M. Groppi, G. Russo, and G. Stracquadanio, High order semi-lagrangian methods for the bgk equation (2014), arXiv:1411.7929 [math.NA] .
  * Krämer  _et al._ [2017] A. Krämer, K. Küllmer, D. Reith, W. Joppich, and H. Foysi, Phys. Rev. E 95 (2017).
  * Krämer  _et al._ [2020] A. Krämer, D. Wilde, K. Küllmer, D. Reith, H. Foysi, and W. Joppich, Comput. Math. with Appl. 79, 34 (2020).
  * Wilde  _et al._ [2020] D. Wilde, A. Krämer, D. Reith, and H. Foysi, Phys. Rev. E 101, 1 (2020), 1910.13918 .
  * Stroud [1973] A. H. Stroud, Math. Comput. 10.2307/2005635 (1973).
  * Cools [2003] R. Cools, Journal of Complexity 19, 445 (2003), oberwolfach Special Issue.
  * Wilde  _et al._ [2021] D. Wilde, A. Krämer, M. Bedrunka, D. Reith, and H. Foysi, Journal of Computational Science 51, 101355 (2021), arXiv:2101.03982 .
  * Krüger  _et al._ [2017] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M. Viggen,  _The lattice Boltzmann method: Principles and practice_ (Springer Nature, 2017).
  * Chen  _et al._ [2014] L. Chen, Q. Kang, Y. Mu, Y. L. He, and W. Q. Tao, Int. J. Heat and Mass Transfer 76, 210 (2014).
  * Wang  _et al._ [2016] L. P. Wang, C. Peng, Z. Guo, and Z. Yu, Comput. Fluids 124, 226 (2016).
  * Peng  _et al._ [2003] Y. Peng, C. Shu, and Y. T. Chew, Phys. Rev. E 68, 8 (2003).
  * Geier  _et al._ [2015] M. Geier, M. Schönherr, A. Pasquali, and M. Krafczyk, Comp. and Math. with Appl. 70, 507 (2015).
  * Dorschner  _et al._ [2016] B. Dorschner, F. Bösch, S. S. Chikatamarla, K. Boulouchos, and I. V. Karlin, J. Fluid Mech. 801, 623 (2016).
  * Alexander  _et al._ [1992] F. J. Alexander, H. Chen, S. Chen, and G. D. Doolen, Phys. Rev. A 46, 1967 (1992).
  * Frapolli  _et al._ [2015] N. Frapolli, S. S. Chikatamarla, and I. V. Karlin, Phys. Rev. E 92, 061301 (2015).
  * Frapolli  _et al._ [2016a] N. Frapolli, S. S. Chikatamarla, and I. V. Karlin, Phys. Rev. Lett. 117, 10604 (2016a).
  * Frapolli  _et al._ [2016b] N. Frapolli, S. S. Chikatamarla, and I. V. Karlin, Phys. Rev. E 93, 063302 (2016b).
  * Feng  _et al._ [2016] Y. Feng, P. Sagaut, and W. Q. Tao, Comput. Fluids 131, 45 (2016).
  * Coreixas  _et al._ [2017] C. Coreixas, G. Wissocq, G. Puigt, J. F. Boussuge, and P. Sagaut, Phys. Rev. E 96 (2017).
  * Atif  _et al._ [2018] M. Atif, M. Namburi, and S. Ansumali, Phys. Rev. E 98, 053311 (2018).
  * Dorschner  _et al._ [2018] B. Dorschner, F. Bösch, and I. V. Karlin, Phys. Rev. Lett. 121, 130602 (2018).
  * Feng  _et al._ [2019] Y. Feng, P. Boivin, J. Jacob, and P. Sagaut, J. Comput. Phys. 394, 82 (2019).
  * Hosseini  _et al._ [2019] S. A. Hosseini, C. Coreixas, N. Darabiha, and D. Thévenin, Phys. Review E 100, 1 (2019).
  * Saadat  _et al._ [2019] M. H. Saadat, F. Bösch, and I. V. Karlin, Phys. Rev. E 99, 013306 (2019).
  * Saadat  _et al._ [2020] M. H. Saadat, F. Bösch, and I. V. Karlin, Phys. Rev. E 101, 23311 (2020).
  * Latt  _et al._ [2020] J. Latt, C. Coreixas, J. Beny, and A. Parmigiani, Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 378, 20190559 (2020).
  * Farag  _et al._ [2020] G. Farag, S. Zhao, T. Coratger, P. Boivin, G. Chiavassa, and P. Sagaut, Phys. Fluids 32, 10.1063/5.0011839 (2020).
  * Coreixas and Latt [2020] C. Coreixas and J. Latt, Physics of Fluids 32, 116102 (2020), https://doi.org/10.1063/5.0027986 .
  * Shan and He [1998] X. Shan and X. He, Phys. Rev. Lett. 80, 65 (1998), 9712001 [comp-gas] .
  * He  _et al._ [1998] X. He, S. Chen, and G. D. Doolen, J. Comput. Phys. 146, 282 (1998).
  * Prasianakis and Karlin [2007] N. I. Prasianakis and I. V. Karlin, Phys. Rev. E 76, 016702 (2007).
  * Bardow  _et al._ [2006] A. Bardow, I. V. Karlin, and A. A. Gusev, Europhys. Lett. 75, 434 (2006).
  * Lee and Lin [2003] T. Lee and C.-L. L. Lin, J. Comput. Phys. 185, 445 (2003).
  * Min and Lee [2011] M. Min and T. Lee, J. Comput. Phys. 230, 245 (2011).
  * Guo  _et al._ [2013] Z. Guo, K. Xu, and R. Wang, Phys. Rev. E. 88, 33305 (2013).
  * Guo  _et al._ [2014] Z. Guo, R. Wang, and K. Xu, Phys. Rev. E 91, 033313 (2014), 1406.5668 .
  * Arasaratnam and Haykin [2009] I. Arasaratnam and S. Haykin, IEEE Trans. on Automatic Control 54, 1254 (2009).
  * Thantanapally  _et al._ [2013] C. Thantanapally, S. Singh, D. V. Patil, S. Succi, and S. Ansumali, International Journal of Modern Physics C 24, 1340004 (2013).
  * Chapman and Cowling [1970] S. Chapman and T. G. Cowling,  _The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases_ (Cambridge University Press, Cambridge, 1970).
  * Wilde  _et al._ [2019] D. Wilde, A. Krämer, K. Küllmer, H. Foysi, and D. Reith, International Journal for Numerical Methods in Fluids 90, 156 (2019).
  * Shan  _et al._ [2006] X. Shan, X. F. Yuan, and H. Chen, J. Fluid Mech. 550, 413 (2006).
  * Haegemans and Piessens [1977] A. Haegemans and R. Piessens, SIAM J. on Num. Analys. 14, 492 (1977), https://doi.org/10.1137/0714029 .
  * Konyaev [1977] S. Konyaev, in  _Doklady Akademii Nauk_ , Vol. 233 (Russ. Acad. of Sci., 1977) pp. 784–787.
  * Van Zandt [2019] J. R. Van Zandt, ETNA - Electron. Trans. Numer. Anal. 51, 219 (2019).
  * Uga  _et al._ [2013] K. C. Uga, M. Min, T. Lee, and P. F. Fischer, Comp. and Math. with Appl. 65, 239 (2013).
  * Lax and Liu [1998] P. Lax and X. Liu, SIAM Journal on Scientific Computing 19, 319 (1998).
  * Kurganov and Tadmor [2002] A. Kurganov and E. Tadmor, Numer. Methods Partial Differ. Equ. 18, 584 (2002).
  * Peng and Yang [2018] N. Peng and Y. Yang, Physical Review Fluids 3, 1 (2018).
  * Lusher and Sandham [2021] D. J. Lusher and N. D. Sandham, AIAA J. 59, 533 (2021).
  * Pruett [2008] C. Pruett, Theor. Comput. Fluid Dyn. 22, 275 (2008).
  * Oberle  _et al._ [2020] D. Oberle, C. D. Pruett, and P. Jenny, Phys. Fluids 32, 10.1063/5.0006637 (2020).
  * Siebert  _et al._ [2008] D. N. Siebert, L. A. Hegele, and P. C. Philippi, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 77, 10.1103/PhysRevE.77.026707 (2008).
  * Sun and Hsu [2003] C. Sun and A. T. Hsu, Phys. Rev. E - Stat. Physics, Plasmas, Fluids, Relat. Interdiscip. Top. 68, 14 (2003).
  * Watari and Tsutahara [2006] M. Watari and M. Tsutahara, Phys. A Stat. Mech. its Appl. 364, 129 (2006).
  * Noah  _et al._ [2021] K. Noah, F. S. Lien, and E. Yee, Int. J. Numer. Methods Fluids 93, 927 (2021).
  * Chen  _et al._ [2020] T. Chen, X. Wen, L. P. Wang, Z. Guo, J. Wang, and S. Chen, Phys. Fluids 32, 10.1063/5.0029424 (2020).
  * Pavlo  _et al._ [1998] P. Pavlo, G. Vahala, and L. Vahala, Phys. Rev. Lett. 80, 3960 (1998).
  * Pavlo  _et al._ [2002] P. Pavlo, G. Vahala, and L. Vahala, J. Stat. Phys. 107, 499 (2002).
  * Zakirov  _et al._ [2019] A. V. Zakirov, B. A. Korneev, V. D. Levchenko, and A. Y. Perepelkina, Keldysh Inst. Prepr. , 1 (2019).
  * Frapolli [2017] N. Frapolli,  _Entropic lattice Boltzmann models for thermal and compressible flows_ , Ph.D. thesis, ETH Zürich, Switzerland (2017).
  * Frapolli  _et al._ [2020] N. Frapolli, S. Chikatamarla, and I. Karlin, Entropy 22, 10.3390/e22030370 (2020).
  * Nie  _et al._ [2009] X. Nie, X. Shan, and H. Chen, 47th AIAA Aerosp. Sci. Meet. Incl. New Horizons Forum Aerosp. Expo. 10.2514/6.2009-139 (2009).
  * Renard  _et al._ [2021] F. Renard, Y. Feng, J.-F. Boussuge, and P. Sagaut, Computers & Fluids 219, 104867 (2021).
  * Saadat  _et al._ [2021] M. H. Saadat, S. A. Hosseini, B. Dorschner, and I. V. Karlin, Phys. Fluids 33, 046104 (2021), arXiv:2102.09512 .
  * Kronbichler and Kormann [2019] M. Kronbichler and K. Kormann, ACM Transactions on Mathematical Software 45, 10.1145/3325864 (2019), arXiv:1711.03590 .
  * Küllmer [2019] K. Küllmer,  _Pseudopotentialbasierte Lattice-Boltzmann-Verfahren für Mehrphasenströmungen_ , Ph.D. thesis, Universität Siegen, Germany (2019).
  * Coreixas [2018] C. Coreixas,  _High-order extension of the recursive regularized lattice Boltzmann method_ , Ph.D. thesis, INP Toulouse, France (2018).