Semi-Lagrangian lattice Boltzmann method for compressible flows

Dominik Wilde ,1,2,* Andreas Krämer ,3 Dirk Reith ,2,4 and Holger Foysi1 1Department of Mechanical Engineering, University of Siegen, Paul-Bonatz-Straße 9-11, D-57076 Siegen-Weidenau, Germany 2Institute of Technology, Resource and Energy-efficient Engineering (TREE), Bonn-Rhein-Sieg University of Applied Sciences, Grantham-Allee 20, D-53757 Sankt Augustin, Germany 3National Heart, Lung, and Blood Institute, National Institutes of Health, Bethesda, Maryland 20892, USA 4Fraunhofer Institute for Algorithms and Scientific Computing (SCAI), Schloss Birlinghoven, D-53754 Sankt Augustin, Germany


I. INTRODUCTION
In the field of weakly compressible and isothermal flows, the lattice Boltzmann method (LBM) [1][2][3] has emerged as an efficient numerical solver that suits modern, highly parallel computing architectures. Consequently, many attempts have been undertaken to extend the method to energy conserving flows [4]. Despite considerable progress in this field, the research on robust and efficient lattice Boltzmann models for compressible flows is still ongoing.
Existing approaches can generally be categorized in two ways [5]: on-versus off-lattice, and single-versus doublepopulation. While on-lattice approaches inherit the LBM's exact, space-filling streaming step [6][7][8][9], off-lattice methods discretize the Boltzmann equation through finite volume [10][11][12][13] or finite difference schemes [14][15][16]. Orthogonal to this classification, single-population models [4,17] express all physical moments, including energy and heat flux, by a single distribution function. In order to represent all moments, the particle velocity sets have to be extended beyond the lattices typically used for weakly compressible flows. In contrast, double-population models [8,18,19] represent the local internal energy through a separate distribution function that is coupled with the density and momentum coming from the first distribution function. * wilde.aerospace@gmail.com The works of Frapolli [17] and Coreixas [20] present an excellent overview of compressible extensions of the LBM. They show how previous approaches have traditionally suffered from instability, numerical dissipation, large computational cost, and undesirable couplings between the physical and numerical parameters.
The use of an efficient off-lattice method is key, since compressible on-lattice Boltzmann methods suffer from decisive disadvantages. First, the time step size is set to unity and cannot be changed independently from the physical parameters. In the weakly compressible regime, the Mach number is merely a numerical parameter following the law of similarity, so that the time step size can be changed via the flow velocity [5]. In compressible flows, however, the Mach number regains its meaning as a physical quantity with only limited remaining options to change the time step, e.g., by increasing or decreasing the spatial resolution, which is still coupled to the time step size. Alternatively, the speed of sound can be varied, at the price of an inferior numerical stability [4]. Second, the simulation of energy conserving flows requires disproportionately large velocity sets since integer-based values of the particle velocities lead to leading-order errors in the relevant moments of the Maxwell-Boltzmann distribution function [21]. A common way to reduce these errors is to increase the number of abscissae, leading to velocity sets like D1Q7, D2Q49, or D3Q343 [22]. For weakly compressible flows, the D1Q3 velocity set and its tensor products D2Q9 and D3Q27 in combination with a second-order truncation of the Maxwell-Boltzmann distribution mostly match the desired properties of the Navier-Stokes equations in the hydrodynamic limit, except the "famous O(Ma 3 ) error term" [23], which breaks Galilean invariance and restricts simulations to low Mach numbers.
The physical and numerical parameters can be decoupled in an off-lattice framework, where the velocity set does not need to match the computational grid. In the field of off-lattice Boltzmann methods the semi-Lagrangian lattice Boltzmann method (SLLBM) was recently introduced [24] and further investigated in a number of subsequent works [25][26][27][28][29]. To formulate the streaming step the algorithm follows the trajectories of the lattice Boltzmann equation along their characteristics to find the cells of the corresponding departure points. In each cell, a finite element formulation reconstructs the distribution function values from the local degrees of freedom by interpolation polynomials. The SLLBM overcomes the aforementioned drawbacks of the standard lattice Boltzmann method. It decouples the time step from the grid spacing and allows for high-order spatial accuracy on irregular grids. In common with Eulerian-based off-lattice Boltzmann solvers the SLLBM enables the use of non-integer-based velocity sets. Yet, compared to Eulerian-based off-lattice Boltzmann schemes, the streaming step is computationally more efficient [28], reduces numerical dissipation, and is still close to the original formulation of the LBM, as it also follows the trajectories of the discrete Boltzmann equation back in time [24]. The problems that arise with interpolation-based schemes have been discussed in earlier works, e.g., by Chen [10], like the issue that interpolation-based LBM do not conserve mass, momentum, and energy. However, this effect can be at least reduced by using high-order interpolation polynomials, which we showed in a recent work [28].
The present work thoroughly investigates an SLLBM solver for compressible flows. The methodology is related to the Particles on Demand (PoD) method [25], which attained considerable interest [29]. The PoD made use of a semi-Lagrangian streaming step to attain the off-lattice distribution function values in a dynamically shifting reference frame, similar to other adaptive LBM models for compressible flows [30][31][32]. In contrast to PoD, we restrict the present method to static, nonmoving reference frames at a fixed reference temperature and show that the method still remains competitive in terms of Galilean invariance up to a certain Mach number. In that sense, the present method is closer to the original lattice Boltzmann method formulation, which also operates in a static, nonmoving reference frame.
In addition to the achievements of PoD, we use Hermitebased equilibria [33,34] and high-order interpolation polynomials. We apply the double-population off-lattice approach that mitigates the aforementioned issues of previous solvers and provides a more comprehensive solution for simulating compressible flows with the LBM.
The semi-Lagrangian streaming step also offers more flexibility in constructing velocity sets from Hermite polynomials [35,36], which is particularly relevant for compressible flow. Originating in Grad's work on the moment system in kinetic theory [37], a number of subsequent works showed that consistent equilibrium distribution functions can be derived by projecting the Boltzmann equation onto a low-dimensional Hilbert subspace through Hermite expansion [21,34,36]. This method is used to construct high-order velocity sets that recover the compressible Fourier-Navier-Stokes equations in the hydrodynamic limit [23,38].
This article is structured as follows. Section II derives the equilibria and velocity sets and recapitulates the semi-Lagrangian advection scheme. Section III validates the method by simulations of a moving incompressible Taylor-Green vortex, a smooth density propagation, a Sod shock tube, a two-dimensional Riemann problem, and a shock-vortex interaction. Section IV discusses the properties of the compressible SLLBM and highlights the advantages over other compressible LBM solvers. Section V provides conclusions.

II. METHODOLOGY
The Boltzmann equation reads with the particle distribution f , the particle velocities ξ , and the collision operator . In this work, only the linearized single relaxation time collision operator proposed by Bhatnagar-Gross-Krook (BGK) [39] will be used: where f eq denotes the equilibrium distribution function, while the relaxation time λ = ν/c 2 s reflects the kinematic viscosity ν and the speed of sound c s . The macroscopic density (ρ), velocity (u), and energy (E ) are represented by the moments of the distribution function: 2ρE = ( f |ξ| 2 + g)dξ.
The second distribution function g represents the rotational energy of the polyatomic molecules and is described in Sec. II D.

A. Equilibrium distribution function
One key component to correctly calculate thermal and compressible flows is the construction of the equilibrium distribution function f eq . The present section is based on the comprehensive derivation of the equilibrium, which was detailed by Shan et al. [34]. The equilibrium is found via scaled Hermite polynomials H (n) (ξ), which are formulated in terms of the normalized particle velocitiesξ = ξ/c s . The Hermite polynomials are defined by the generating function, as where D is the dimension and c 2 s = RT 0 is the product of gas constant and reference temperature. (Note that the generating function is differentiated with respect to ξ rather thanξ.)  The scaled Hermite polynomials up to fourth order are Here, I is the identity matrix, ξ n := ξ ⊗ ξ ⊗ · · · ⊗ ξ (n times) and for arbitrary vectors, for example, we write ab := a ⊗ b. Furthermore, we defined for an nth order tensor T (n) , according to [40], the relation, with π indicating the significant permutations of either the base vectors or indices. For arbitrary vectors v andξ one has, for example, vξ 2 = 1 3 (vξξ +ξvξ +ξξv). The (equilibrium) distribution functions can be expressed as a Hermite series [23] of order N, with : indicating full contraction. The order N determines the physics of the lattice Boltzmann model. The Hermite coefficients a are obtained via projection of f on the orthogonal basis of Hermite tensors. The equilibrium coefficients a (n) eq are therefore directly related to the moments of the Maxwell-Boltzmann equilibrium, as [23] a (0) eq = ρ, (10a) a (4) αβγ δ,eq = R eq αβγ δ = ρ u α u β u γ u δ + T 0 (θ − 1)((δ αβ δ γ δ +δ αγ δ βδ + δ αδ δ βγ )T 0 (θ − 1) Note that in case of the established weakly compressible lattice Boltzmann simulations, the isothermal assumption leads to a local temperature of θ = T /T 0 = 1 and to vanishing terms in Eq. (10). Truncating at second order covers the well-known lattice Boltzmann model for weakly compressible flows at low velocity. Truncating at third order compromises the nonequilibrium parts of Q αβγ due to deviations in R eq αβγ δ . Thus, to represent the heat flux correctly, a fourth-order expansion is required [23,34]. Consequently, to recover the full Fourier-Navier-Stokes equations, the approximation in Eqs. (8) and (9) is of fourth order in the present work.

B. Hermite polynomial velocity sets
To recover the Mth order moment a M , the velocity space is discretized by a (weighted) Gauss-Hermite quadrature, provided N M and with the discrete equilibrium, Here, w i represents the weights coming from a Gauss-Hermite quadrature, and V is the number of velocities in the velocity set. In general, the number of possible velocity sets is several times larger than in on-lattice Boltzmann methods, since the abscissae of the velocity sets are not required to match the grid points. The present work uses a D2Q25 velocity set, which is derived by a tensor multiplication from the one-dimensional D1Q5 velocity set [36]. The one-dimensional abscissae and weights are listed in Table I, while the two-dimensional D2Q25 velocity set is depicted in Fig. 1. Unless scaled, the characteristic lattice speed of sound equals c s = 1. This off-lattice velocity set is of quadrature order Q = 9 and consequently recovers moments up to order N = 4, fulfilling the requirement (M + N )/2 Q to recover the Fourier-Navier-Stokes equations [41].

C. Variable Prandtl number
In the context of thermal and compressible flows, using the BGK operator limits the Prandtl number Pr, i.e., the ratio of kinematic viscosity ν and thermal diffusivity α, to Pr = ν/α = 1. This drawback can be circumvented by using a multirelaxation time collision operator to decouple both properties of the fluid [42], a quasiequilibrium [6,43]  using a BGK-Shakhov model [44]. Although we observed that for the test cases in this work there is no significant difference, the SLLBM solver is capable to adapt the Prandtl number by the quasiequilibrium f * i and g * i [6,43], which is with the nonequilibrium part of the centered heat flux tensor with the same proceeding for g * i .

D. Variable heat capacity ratio
As the Boltzmann equation and the Maxwell-Boltzmann distribution are valid in the assumption of monoatomic flows, the heat capacity ratio γ = C p /C v is fixed to γ = (D + 2)/2, with the dimension D and the heat capacities at constant pressure and volume C p and C v , respectively. Nonetheless, several attempts in the literature solve this issue by using a second set of distribution functions [17,45], inspired by the Rykov model in kinetic theory for polyatomic gases [46], or by adapting the equilibrium distribution to the heat capacity of the fluid [47]. We adopt the ansatz by Frapolli et al. for the LBM by incorporating an additional set of distributions g and its equilibrium g eq .
The equilibrium of g eq can easily be determined as with γ = R/C v + 1, where the gas constant R is usually set to unity (without loss of generality). The distribution function g i follows the same stream and collide algorithm as the distribution function f i and is ultimately used in Eq. (5).
The lattice Boltzmann equations of the present model read with the relaxation time τ = ν/(c 2 s δ t ) + 0.5, where the term of 0.5 is a consequence of the second-order time integration in terms of the trapezoidal rule [48][49][50]. The relaxation time τ Pr = (τ − 0.5)/Pr + 0.5 is related to the Prandtl number Pr. In contrast to the usual lattice Boltzmann algorithm, neither the time step size δ t , nor the expression δ t ξ i are required to be integers due to the use of a semi-Lagrangian streaming step.

E. Semi-Lagrangian streaming step
To obtain the off-vertex distribution function values f i (and likewise g i ), the SLLBM uses a cell-wise polynomial interpolation. Therefore, the computational domain is divided into a finite number of cells N with Lagrange polynomials ψ of order p as shape functions. The distribution function values at the support points in cell are denoted asf i j , so that for all points x in cell .
The distributionsf i j are updated via Eq. (16), where the value at the departure point f i (x − δ t ξ i ) is obtained via the interpolation defined in Eq. (17). An equidistant interpolation is stable for low-order interpolation polynomials only. Therefore, at higher interpolation orders this work applies Gauß-Lobatto-Chebyshev support points to dampen the Runge phenomenon and to increase the stability [24,51]. Gauß-Lobatto-Chebyshev support points in one dimension read The two-dimensional unit cell used in this work is depicted in Fig. 2.
The interpolation error of semi-Lagrangian advection schemes has been investigated, for instance, by Einkemmer and Ostermann [52] or Falcone and Ferretti [53]. For the SLLBM used in the present work, we have previously discussed the nontrivial behavior of the temporal discretization error [28], which is of order, This means in practice that the error has a minimum for intermediate CFL numbers and is bounded by the spatial discretization error in the limit of small time steps.
All SLLBM simulations in this manuscript have been performed with the NATRIUM library [28], which makes use of the finite element library DEAL.II [54]. See [28] for further details concerning the implementation.

F. Thermo-hydrodynamic limit
The thermo-hydrodynamic limit of the presented equilibrium distribution function was presented and thoroughly discussed by Malaspinas et al. [38] and by Coreixas et al. [23], who found with the stress tensor, with the kinematic viscosity, and with the heat conductivity, The bulk viscosity equals [22] with C v and C p being the heat capacity at constant volume and pressure, respectively.

G. On-lattice Boltzmann solver for comparison
Section III C compares the SLLBM with an on-lattice Boltzmann solver. For the latter, we chose the D2V37 velocity set proposed by Philippi et al. [55], which is of the same quadrature order Q = 9 as the D2Q25 in this work. Additionally, the equilibrium in Sec. II A and the second distribution function in Sec. II D were used to obtain a variable adiabatic exponent γ . The on-lattice simulations were performed by the solver LETTUCE [56], which is based on the machine learning framework PYTORCH [57]. The code is written in PYTHON and makes use of the PYTORCH routines to enable GPU simulations.

III. RESULTS
In this section, a variety of test cases will be studied. First, the errors of the method are quantified by simulations of the two-dimensional incompressible Taylor-Green vortex. Then the accuracy will be confirmed by a smooth density propagation. The Sod shock tube shows the general capability of the SLLBM to deal with shocks at the original density ratio presented by Sod [58]. Then the two-dimensional Riemann problem will be studied, followed by a simulation of a shockvortex interaction on nonuniform grids.

A. Moving two-dimensional Taylor-Green vortex
The SLLBM was validated by a moving incompressible two-dimensional Taylor-Green vortex. As it comes along with a reference solution it is suited to quantify the interpolationcaused errors and to show the high degree of Galilean invariance if an appropriate velocity set is used. The initial and reference solution is (2y) We compare two velocity sets: the standard D2Q9 velocity set and the D2Q25 velocity set, both with the SLLBM advection step. The fourth-order equilibrium specified in Sec. II A was used in an isothermal configuration, i.e., θ = 1.    Figure 3(a) shows the nonmoving case. As detailed in Sec. II E the interpolation provoked nonnegligible errors, which were in total two to three times larger for the D2Q25 velocity set in comparison to the D2Q9, caused by the different amount of interpolations needed to perform one time step. However, the larger velocity set pays off for constantly moving vortices in case Fig. 3(b). Here the D2Q9 simulation converged towards an error, which is caused by the lacking Galilean invariance of this particular velocity set. In contrast, the D2Q25's errors were nearly unchanged compared to simulations [ Fig. 3

(a)] with u h = 0.
This test case also served to evaluate the maximum Mach number ranges of the D2Q25 velocity set. To achieve this, the horizontal velocity Ma h was increased until the simulations became unstable, while the initial vortex velocity remained at Ma v = 0.01. Figure 4 depicts the numerical errors over the horizontal Mach number Ma h for three different time step sizes with two main findings. First, there is a local minimum of the numerical error in case of the time step size δ t = 0.0004. Second, for all three time step sizes the numerical errors remain more or less constant over the range of horizontal Mach numbers, while the Mach number limit was abruptly exceeded. Obviously, the maximum Mach number of semi-Lagrangian schemes is not solely determined by equilibrium and velocity set, but also by the time step size.

B. Smooth density propagation
To measure the accuracy of the SLLBM for compressible flows, the smooth density propagation [59,60] was evaluated.
The exact solution is The domain was x, y ∈ [−1, 1] and the simulation end time was t = 2. For the coarsest grid of 16 × 16 grid points (i.e., 4 × 4 cells and p = 4), this time was reached in 1105 steps; the number of steps doubled for each finer simulation. The kinematic viscosity was set to ν = 10 −10 ; the Mach number was set to Ma = 0.2. Table II depicts the numerical errors of the density for six different resolutions. At coarse resolutions the numerical errors decreased with approximately fourth order. At the finest resolutions, the second-order temporal order became dominant.

C. Sod shock tube
In a first test case with discontinuities, the well-known Sod shock tube [58] was simulated. The one-dimensional domain x ∈ [0, 1] was divided into two regions at x = 0.5 with the initial conditions,  Figure 5 depicts the density, velocity, and temperature for the polynomial order p = 4 at t = 0.15, which is reached after 735 time steps.
The results accurately matched the reference solution [58] with only minor oscillations visible, comparable to other recently introduced methods [60].
The proposed method in our tests demonstrated clearly improved numerical stability compared to previous on-lattice methods. To compare the stability properties of the SLLBM and of the representative on-lattice Boltzmann solver, specified in Sec. II G, the Sod shock tube was again simulated for a variety of resolutions and physical viscosities, ν phys , until t = 1.5. A simulation was rated as stable, if the simulation finished successfully and all macroscopic values were still bounded. Figure 6 shows that the on-lattice Boltzmann solver was unstable for most of the lower viscosities at low resolutions. A simulation at ν = 0.0001, an example of an often used setting for Sod shock tubes in the literature, was only stable for the on-lattice LBM using an excessive resolution of N = 6400. Contrary to that, the flexible time step size of the SLLBM allowed for stable simulations even at low resolutions. Figure 7 depicts all stable simulations for three different time step sizes δ t = 0.25δ t , δ t = 0.50δ t , and δ t = 1.00δ t , where δ t is the time step size of the respective on-lattice simulation. All of the representative resolution and viscosity scenarios were simulated stably with at least one time step size by the SLLBM.   [61].
Taken together, the Sod shock tube shows that the proposed SLLBM solver is capable to accurately capture shocks and convinces by a larger stability region than the D2V37 onlattice Boltzmann simulation.

D. 2D Riemann
The two-dimensional Riemann problem is another demanding test case for compressible inviscid flow solvers, which was deeply investigated for a variety of initial configurations by Lax and Liu [61] and Kurganov and Tadmor [62]. The original works made use of Riemann solvers and semidiscrete finite difference schemes, respectively, to obtain the reference solutions.
The domain (x, y) ∈ [0, 1] was divided into four quadrants with the initial conditions shown in Table III. The simulation domain was discretized into N = 128 × 128 cells with an order of finite elements of p = 4, which results in N = 512 × 512 grid points. The relaxation time τ = 1.94815 was determined according to τ = ν/(c 2 s δ t ) + 0.5. The physical simulation end time was t end = 0.25, which was reached after 2091 iterations, i.e., the viscosity is ν = 0.000173, assuming the reference temperature T 0 = 1. The heat capacity ratio was set to γ = 1.4. At the boundaries, a zero-gradient condition, i.e., ∂ f = 0, was applied. From a practical point of view, to free the simulation domain from disturbing boundary interactions, the domain was set to twice the size and clipped during post-processing. At first, Gauss-Lobatto-Chebyshev nodes were used to distribute the support points within each cell. Figure 8 shows the density contours after 2091 iterations. In total, 45 equally distributed contours were applied. The results were in good agreement with the reference solution in [61,62], although minor oscillations were still visible near the shock fronts. Figure 9 depicts the pressure, velocity, and temperature charts along the y = x diagonal at the simulation end time. Although there is no reference solution available, the plot gives a good impression of the strong discontinuities that are also the main source of instabilities in this test case.
To emphasize the value of nonequidistant support points for the interpolation procedure, the same configuration was simulated again with equidistant nodes. The simulation became unstable shortly after 200 time steps. Figure 10 shows the pressure along the y = x axis and reveals the heavy oscillations of the simulation with equidistant nodes shortly before the simulation diverges.  [61] and [62].

E. Shock vortex interaction
The final test case was a 2D shock vortex interaction, which was originally presented by Inoue and Hattori [63]. For this simulation a transformed grid was used, which shows the applicability of the SLLBM to nonuniform grids. The simulation setup was as follows: The domain consists of two regions, which are separated by a steady shock. The Mach number of region A is Ma A = −1.2, while the Mach number in region B Ma B ≈ −0.842 is determined by the Rankine-Hugoniot condition. A vortex with characteristic Mach number Ma v = 0.5 is advected in region A towards the shock (from right to left) and interacts with it throughout the simulation. The flow field of the vortex is given by The initial pressure and density field are and The temperature field can be determined via the ideal gas law P = ρRT.
The Reynolds number is defined by with R being the resolution of the characteristic radius r. The Prandtl number is Pr = 0.75. The simulation was run on a physical domain of size 60 × 24, which was discretized by N = 256 × 256 cells and the polynomial order was p = 4, resulting in 1024 lattice nodes in each direction. However, to highly resolve the steady shock the grid was transformed in the x direction according to a simple transformation function: where is a stretching parameter, which was set to = 0.95. A similar refinement was done by Inoue and Hattori in the original work [63]. Although the shock slightly moved during the interaction with the shock, the vortex was still located in the refinement zone, which led to satisfactory results and reduced the amount of needed grid points. The shock was located at x shock = 30 and the initial vortex center position was x vortex = 32. We observed that the DNS results by Inoue and Hattori were only matched if the vortex was initialized exactly at this position although partly located in the shock. The schematic grid of the computation is shown in Fig. 11. The average resolution of the characteristic radius along the x axis r was 1024/60 ≈ 17.07 grid points and 1024/24 ≈ 42.67 grid points along the y axis. The large domain and periodic boundary conditions in all directions were sufficient to keep the domain of interest free from disturbing influences of the boundaries. The characteristic time is defined as t = R/c s . The time step size was δ t = 7 × 10 −4 so that the final result at t = 8 was reached in 11 290 iterations. Figure 12 shows the density contours in the range ρ ∈ [0.92, 1.55] with 119 contours. The contours were in good agreement with previous works, e.g., [17,63]. Both the static  shock and the attached shock were well resolved and the density contours appeared without any oscillations.
To quantify the results, the radial sound pressure in a simulation with Ma v = 0.25 and Re = 800 was evaluated and compared with DNS results from Inoue and Hattori [63]. The total mass in the previously presented well-resolved test cases was found to be more or less constant over the entire simulation time. However, since the SLLBM is not rigorously mass-conserving, Figure 14 presents the evolution of mass over 1500 simulation steps for a slightly under-resolved configuration of Re = 800, Ma = 0.25, and 512 × 512 grid points at two polynomial orders. In case of a low-order polynomial interpolation of p = 2, there was a minor but steady loss of mass, which amounted to approximately 0.05% of the initial mass after 1500 steps. The mass conservation was tighter for the high-order interpolation. For p = 4, the mass oscillated around the initial value as opposed to a systematic gain or loss. 053306-9

IV. DISCUSSION
As detailed in the introduction of this work, the research on lattice Boltmann methods for compressible flows is ongoing. We have shown how the use of an efficient off-lattice streaming step enables compressible simulations by relaxing the coupling between time step, computational grid, and particle velocities in the LBM.
Originating from the usual lattice Boltzmann equation the SLLBM algorithm follows the trajectories along the characteristics and recovers the off-lattice departure points using a finite element interpolation. In common with finite volume or finite difference LBMs, this procedure allows for a flexible time step size. This key feature is missing in usual on-lattice Boltzmann solvers for compressible flows. Contrary to the Eulerian methods, the SLLBM avoids costly and often highdissipative time integration. The compressible SLLBM is related to the PoD method, but operates in a static, nonmoving reference frame. Despite this regression, the SLLBM was capable to simulate a variety of compressible test cases. As the weights in the present model are not temperature-dependent, the temperature ratio is-in principle-not limited. Simple temperature diffusion tests not shown in this article revealed a stable temperature ratio of at least 1:1000.
We tested the compressible SLLBM for five common test cases, in both inviscid and viscous regimes. The twodimensional Taylor-Green vortex showed the high degree of Galilean invariance of the scheme. The accuracy test by a smooth density propagation confirmed the fourth order in space and second order in time of the scheme. The Sod shock tube in its original configuration and a pressure ratio of P 1 /P 2 = 10 were successfully tested with the SLLBM. In contrast, on-lattice Boltzmann solvers like the one presented in [8] often restrict the pressure ratio to lower values on larger grids. In contrast, the SLLBM's time step size can be tuned to the needs of the simulation, even with the most common BGK model. Furthermore, we employed fourth-order interpolation polynomials, rendering the solution in the smooth areas spatially to the corresponding order [24], which is another unique property of this scheme. However, as shown in the twodimensional Riemann problem in Sec. III D, leveraging the spatial order is only workable for an appropriate distribution of support points within each cell. While equidistant support points led to unstable simulations, the use of fourth-order Gauss-Lobatto-Chebyshev points enabled stable simulations. Moreover, the simulation of the shock-vortex interaction showed that a high-order interpolation improves mass conservation, which is not rigorously enforced in interpolationbased LBM [29]. The equilibrium distribution function based on a Hermite expansion of fourth order is computationally costlier than the usual second order expansion, but conserves energy as required for compressible flows [34]. Concerning the discretization of the velocity space, a D2Q25 velocity set based on the fifth-order Hermite polynomial was used, which is suited to recover the moments up to fourth order correctly [21]. In contrast, on-lattice velocity sets usually provoke deviations in the recovery of the high-order moments [21]. The combination of the equilibrium and velocity set allowed for a total Mach number of Ma = 1.7 in the simulation of the shock-vortex interaction, but the results of the 2D Taylor-Green vortex flow suggest that even higher Mach numbers are attainable, depending on the time step size. A stability analysis [64,65] could clarify the exact limits of the scheme. An obvious way to further increase the Mach number is to incorporate statically shifted [22,66] or dynamically shifted lattices [25]. Future works should thoroughly investigate the path to extend the velocity sets to higher velocities in three dimensions, while keeping the number of discrete velocities lower than 125. The effects of lattice pruning, which were thoroughly studied in the weakly compressible case, e.g., D3Q27 to D3Q19 [67], should also be examined for these high-order velocity sets. Additionally, we propose to pursue an approach that was originally supposed by Watari and Hattori [68] for finite difference LBMs. The authors combined platonic solids, dodecahedra, and icosahedra, to construct highly isotropic velocity sets at different velocity shells. It seems to be worth investigating if these velocity sets allow increasing the Mach number in SLLBM simulations. Corresponding results will be published in an upcoming article.
Another issue is the second distribution function value, needed to correctly represent the heat capacity ratio γ . This model was also successfully used by other authors [6,8,17]. A reduction of the velocity set size of g could possibly make the compressible LBM more attractive.
Taken together, the presented SLLBM opens an alternative and attractive field of solvers for compressible flows and, at the same time, combines many insights from previous LBM solvers for compressible flows.

V. CONCLUSIONS
In this article we showed that the semi-Lagrangian streaming step enables efficient lattice Boltzmann simulations of compressible flows even in a static, nonmoving reference frame. In contrast to other methods, the vertices are organized in cells, and interpolation polynomials up to fourth order are used to attain the off-vertex distribution function values.
The presented SLLBM solver for compressible flows convinces by six main advantages.
(1) Flexibility of the time step size.
(2) Absence of a dedicated time integrator as the integration is performed along characteristics.
(3) Applicability to stretched or unstructured meshes. (4) Possibility to employ advantageous velocity sets. (5) Limitation of mass loss by high-order polynomials. (6) Fourth order of convergence in space and second order in time.
The test cases of a 2D Taylor Green vortex, a smooth density propagation, a Sod shock tube, a 2D Riemann problem, and a shock-vortex interaction on transformed grids show the general applicability to flows with discontinuities. Thereby, Gauss-Lobatto-Chebyshev support points reduce oscillations near the shock fronts.
In our future work, the presented framework will be applied to three-dimensional test cases with a focus on further increasing the time step size and on reducing the number of discrete velocities. 053306-10