---
abstract: 'We present high-order, fully explicit projective integration schemes for nonlinear collisional kinetic equations such as the BGK and Boltzmann equation. The methods first take a few small (inner) steps with a simple, explicit method (such as direct forward Euler) to damp out the stiff components of the solution. Then, the time derivative is estimated and used in an (outer) Runge-Kutta method of arbitrary order. The procedure can be recursively repeated on a hierarchy of projective levels to construct telescopic projective integration methods. Based on the spectrum of the linearized collision operator, we deduce that the computational cost of the method is essentially independent of the stiffness of the problem: with an appropriate choice of inner step size, the time step restriction on the outer time step, as well as the number of inner time steps, is independent of the stiffness of the (collisional) source term. In some cases, the number of levels in the telescopic hierarchy depends logarithmically on the stiffness. We illustrate the method with numerical results in one and two spatial dimensions.'
author:
- 'Ward Melis [^1]'
- 'Thomas Rey [^2]'
- 'Giovanni Samaey [^3]'
bibliography:
- 'refs.bib'
title: Projective and telescopic projective integration for the nonlinear BGK and Boltzmann equations
---

**Keywords:** `Boltzmann equation, BGK equation, Projective Integration, spectral theory, fast spectral scheme.`

**2010 MSC.** 82B40, 76P05, 65M70, 65M08, 65M12.

Introduction \[sec:introduction\]
=================================

Kinetic equations represent a gas as a set of particles undergoing instantaneous collisions interspersed with ballistic motion [@CIP]. Nowadays, these models appear in a variety of sciences and applications, such as astrophysics, aerospace and nuclear engineering, semiconductors, fusion processes in plasmas, as well as biology, finance and social sciences. The common structure of such equations consists in a combination of a linear transport term with one or more interaction terms, which together dictate the time evolution of the distribution of particles in the (six-dimensional) position-velocity phase space. From a numerical point of view, it is clear that this results in a real challenge, since the computational cost immediately becomes prohibitive for realistic problems [@DimarcoPareschi2014]. Aside from the curse of dimensionality, there are many other difficulties which are specific to kinetic equations. We recall two among the most important ones. The first is the computational cost related to the evaluation of the collision operator, which implies the computation of multidimensional integrals in each point of the physical space [@FiMoPa2006; @PaRuSINUM2000]. The second challenge is represented by the presence of multiple time scales in the collision dynamics, leading to a very small mean free path, at least in parts of the spatial domain. Usually, computational problems exhibit multiple regimes in different regions in space. This requires the development of adapted numerical schemes to avoid the resolution of the stiff dynamics [@DimarcoPareschi2013; @Jin_review; @Jin1999; @BLM2008; @degondrev].

Historically, two different approaches are generally used to tackle kinetic equations numerically: deterministic methods, such as finite volume, semi-Lagrangian and spectral schemes [@DimarcoPareschi2014], and probabilistic methods, such as Direct Simulation Monte Carlo (DSMC) schemes [@bird; @Caflisch98]. Both methodologies have strengths and weaknesses. Deterministic methods can normally reach high orders of accuracy. Nevertheless, stochastic methods are often faster, especially for solving steady problems, but, typically, exhibit lower convergence rates and difficulties in describing non-stationary and slow motion flows. In this paper, we will consider deterministic methods, in which we evaluate the collision operator using a fast spectral method, in the spirit of [@MoPa:2006]. For a comprehensive overview of numerical schemes for collisional kinetic equations, such as equation , we refer to [@DimarcoPareschi2014] and references therein.

In this paper, we are specifically interested in the time discretization of kinetic equations with stiffness arising from multiple time scales in the collision operator. The stiffness is usually characterized by the (small) mean free path ${\varepsilon}$, and becomes infinite when ${\varepsilon}$ tends to zero. In that limit, a limiting macroscopic equation emerges in terms of a few moments of the particle distribution (density, momentum, energy); the full particle distribution then relaxes infinitely quickly to a Maxwellian distribution defined by these low-order moments. There is currently a large research effort in the design of algorithms that are uniformly stable in ${\varepsilon}$ and approach a scheme for the limiting equation when ${\varepsilon}$ tends to 0; such schemes are called asymptotic-preserving in the sense of Jin [@Jin1999]. Again, we refer to the recent review [@DimarcoPareschi2014] for a clear survey on numerical methods for kinetic equations. Here, we briefly review some achievements using different strategies. In [@Jin1999; @Jin2000a], separating the distribution function $f$ into its odd and even parts in the velocity variable results in a coupled system of transport equations where the stiffness appears only in the source term, allowing to use a time-splitting technique with implicit treatment of the source term; see also related work in [@Jin1999; @Klar1999; @Klar1999a]. Implicit-explicit (IMEX) schemes are an extensively studied technique to tackle this kind of problems, see [@ascher1995; @FilbetJin2010] and references therein. Recent results in this setting were obtained by Dimarco et al. to deal with nonlinear collision kernels [@DimarcoPareschi2013], and an extension to hyperbolic systems in a diffusive limit is given in [@Boscarino2013]. A different approach, based on well-balanced methods, was introduced by Gosse and Toscani [@Gosse2003; @Gosse2004], see also [@Buet2007]. When the collision operator allows for an explicit computation, an explicit scheme can be obtained subject to a classical diffusion CFL condition by splitting the particle distribution into its mean value and the first-order fluctuations in a Chapman-Enskog expansion form [@Godillon-Lafitte2005]. Also closure by moments, e.g. [@Coulombel2005], can lead to reduced systems for which time-splitting provides new classes of schemes [@Carrillo2008]. Alternatively, a micro-macro decomposition based on a Chapman-Enskog expansion has been proposed [@Lemou2008], leading to a system of transport equations that allows to design a semi-implicit scheme without time splitting. A non-local procedure based on the quadrature of kernels obtained through pseudo-differential calculus was proposed in [@Besse2010].

A robust and fully explicit method, which allows for time integration of (two-scale) stiff systems with arbitrary order of accuracy in time, is projective integration. Projective integration was proposed in [@Gear2003projective] for stiff systems of ordinary differential equations with a clear gap in their eigenvalue spectrum. In such stiff problems, the fast modes, corresponding to the Jacobian eigenvalues with large negative real parts, decay quickly, whereas the slow modes correspond to eigenvalues of smaller magnitude and are the solution components of practical interest. Projective integration allows a stable yet explicit integration of such problems by first taking a few small (inner) steps using a step size ${\delta t}$ with a simple, explicit method, until the transients corresponding to the fast modes have died out, and subsequently projecting (extrapolating) the solution forward in time over a large (outer) time step of size ${{\Delta t}> {\delta t}}$. In [@Lafitte2012], projective integration was analyzed for kinetic equations with a diffusive scaling. An arbitrary order version, based on Runge-Kutta methods, has been proposed recently in [@LafitteLejonSamaey2015], where it was also analyzed for kinetic equations with an advection-diffusion limit. In [@LafitteMelisSamaey2017], the scheme was used to construct a explicit, flexible, arbitrary order method for general nonlinear hyperbolic conservation laws, based on relaxation to a kinetic equation. Alternative approaches to obtain a higher-order projective integration scheme have been proposed in [@Lee2007; @Rico-Martinez]. These methods fit within recent research efforts on numerical methods for multiscale simulation [@E2007; @Kevrekidis2003].

For problems exhibiting more than a single fast time scale, telescopic projective integration (TPI) was proposed [@Gear2003telescopic]. In these methods, the projective integration idea is applied recursively. Starting from an inner integrator at the fastest time scale, a projective integration method is constructed with a time step that corresponds to the second-fastest time scale. This projective integration method is then considered as the inner integrator of a projective integration method on yet a coarser level. By repeating this idea, TPI methods construct a hierarchy of projective levels in which each outer integrator step on a certain level serves as an inner integrator step one level higher. The idea was studied and tested for linear kinetic equations in [@MelisSamaey2017]. These methods turn out to have a computational cost that is essentially independent of the stiffness of the collision operator.

We do not call projective integration methods asymptotic-preserving as such, because we cannot explicitly evaluate the scheme for ${\varepsilon}=0$ to obtain a classical numerical scheme for the limiting equation. Nevertheless, projective and telescopic projective integration methods share important features with asymptotic-preserving methods. In particular, their computational cost does (in many cases) not depend on the stiffness of the problem. To be specific, it was shown in [@MelisSamaey2017], for linear kinetic equations, that the number of inner time steps at each level of the telescopic hierarchy is independent of the small-scale parameter ${\varepsilon}$, as is the step size of the outermost integrator. The only parameter in the method that may depend on $epsi$ is the *number* of levels in the telescopic hierarchy. For systems in which the spectrum of the collision operator fall apart into a set of clearly separated clusters (each corresponding to a specific time scale), the number of levels equals the number of spectral clusters. In this situation, the computational cost is completely independent of ${\varepsilon}$. When the collision operator represents a continuum of time scales, the number of projective integration levels increases logarithmically with ${\varepsilon}$.

In this paper, we construct and evaluate telescopic projective integration methods for nonlinear Boltzmann BGK and Boltzmann kinetic equations. The methods are of arbitrary order in time, fully explicit, and general (they do not exploit any particular form of the collision operator). The remainder of this paper is structured as follows. In Section \[sec:models\], we start by presenting the Boltzmann and BGK equations that will be the subject of our simulations. We describe the different projective and telescopic projective integration methods in detail in Section \[sec:time\_integrator\]. (The spatial and velocity discretizations are standard. To make the manuscript self-contained, we present the corresponding numerical methods in Appendices \[sub:weno\] and \[sub:FastSpectral\].) We discuss in Section \[sec:linearOpSpecProp\] the spectral properties of the linearized collision operators, which will guide the choice of the method parameters, ensuring stability of the time integrators. Some numerical experiments are done in Section \[sec:results\] to verify the theory developed in the two previous sections. We conclude in Section \[sec:conclusions\].

Model equations {#sec:models}
===============

In this article, we are interested in rarefied, collisional gases, and then we shall consider Boltzmann-like, collisional kinetic equations. We refer the reader to the classical works [@CIP; @Villani:2002handbook] and the references therein for a more detailed introduction on this vast topic. For a given non-negative initial condition $f_0$, we will study a particle distribution function ${f^{\varepsilon}}= {f^{\varepsilon}}({\mathbf{x}},{\mathbf{v}}, t)$, for $t \geq 0$, ${\mathbf{x}}\in \Omega \subset \mathbb{R}^{D_x}$ and ${\mathbf{v}}\in \mathbb{R}^{D_v}$, solution to the initial-boundary value problem $$\label{eqCollision} 
    {\mathopen{}\mathclose\bgroup\originalleft}\{
    \begin{aligned}
        & \frac{\partial {f^{\varepsilon}}}{\partial t} + {\mathbf{v}}\cdot \nabla_{\mathbf{x}}{f^{\varepsilon}}= \frac{1}{{\varepsilon}} {\mathcal{Q}}({f^{\varepsilon}}), \\
        &\, \\
        & {f^{\varepsilon}}({\mathbf{x}}, {\mathbf{v}}, 0) = f_{0}({\mathbf{x}},{\mathbf{v}}).
    \end{aligned}
    {\aftergroup\egroup\originalright}.$$ The left hand side of equation corresponds to a linear transport operator that comprises the convection of particles in space, whereas the right hand side contains the collision operator that entails velocity changes due to particle collisions. We postpone the description of boundary conditions until Section \[sec:results\], where we discuss the numerical results. We assume that the collision operator fulfils the following three assumptions:

1.  \[hypConservations\] Conservation of mass, momentum and kinetic energy: $$\int_{{\mathbb{R}}^{D_v}} {\mathcal{Q}}(f)({\mathbf{v}}) \, d{\mathbf{v}}= 0, \quad  \int_{{\mathbb{R}}^{D_v}} {\mathcal{Q}}(f)({\mathbf{v}}) \, {\mathbf{v}}\, d{\mathbf{v}}= \bm{0}_{{\mathbb{R}}^{D_v}}, \quad \int_{{\mathbb{R}}^{D_v}} {\mathcal{Q}}(f)({\mathbf{v}}) \, |{\mathbf{v}}|^2 \, d{\mathbf{v}}= 0;$$

2.  \[hypEntropy\] Dissipation of the Boltzmann entropy (H-theorem): $$\int_{{\mathbb{R}}^{D_v}} {\mathcal{Q}}(f)({\mathbf{v}}) \, \log(f)({\mathbf{v}}) \, d{\mathbf{v}}\, \leq \, 0;$$

3.  \[hypEquilib\] Its equilibria are given by Maxwellian distributions: $${\mathcal{Q}}(f) \, = \, 0 \quad \Leftrightarrow \quad f = {\mathcal{M}}_{\mathbf{v}}^{\rho, {\mathbf{{\bar{v}}}}, T} := \frac{\rho}{(2 \pi T)^{D_v/2}} \exp {\mathopen{}\mathclose\bgroup\originalleft}( - \frac{|{\mathbf{v}}-{\mathbf{{\bar{v}}}}|^2}{2 T} {\aftergroup\egroup\originalright}),$$ where the *density* $\rho$, *velocity* ${\mathbf{{\bar{v}}}}$ and *temperature* $T$ of the gas are computed from the distribution function $f$ as: $$\label{eq:f_moments} 
            \rho = \int_{{\mathbb{R}}^{D_v}} f({\mathbf{v}})\,d{\mathbf{v}}, \quad {\mathbf{{\bar{v}}}}= \frac{1}{\rho}\int_{{\mathbb{R}}^{D_v}} {\mathbf{v}}f({\mathbf{v}}) \, d{\mathbf{v}}, 
            \quad T = \frac{1}{D_v \rho} \int_{{\mathbb{R}}^{D_v}} \vert {\mathbf{{\bar{v}}}}- {\mathbf{v}}\vert^2 f({\mathbf{v}}) \,d{\mathbf{v}}.$$

Equation with assumptions \[hypConservations\]-\[hypEntropy\]-\[hypEquilib\] describes numerous models such as the Boltzmann equation for elastic collisions [@Villani:2002handbook] or Fokker-Planck-Landau type equations[@alexandre2000entropy]. The parameter ${\varepsilon}> 0$ is the dimensionless Knudsen number, that is, the ratio between the mean free path of particles and the length scale of observation. It determines the regime of the gas flow, for which we roughly identify the hydrodynamic regime $({\varepsilon}\le 10^{-4})$, the transitional regime $({\varepsilon}\in [10^{-4},10^{-1}])$, and the kinetic regime $({\varepsilon}\ge 10^{-1})$. Moreover, according to assumptions \[hypEntropy\]-\[hypEquilib\], when ${\varepsilon}\to 0$, the distribution ${f^{\varepsilon}}$ converges (at least formally) to a Maxwellian distribution, whose moments are solution to the compressible Euler system for perfect gases, given by: $$\label{eqHydroClosedEuler}
    {\mathopen{}\mathclose\bgroup\originalleft}\{
    \begin{aligned}
        & \partial_t \rho  + \operatorname{div}_{\mathbf{x}}(\rho \, {\mathbf{{\bar{v}}}}) = 0,                \\
        &\, \\
        & \partial_t(\rho \, {\mathbf{{\bar{v}}}}) + \operatorname{div}_{\mathbf{x}}{\mathopen{}\mathclose\bgroup\originalleft}(\rho \, {\mathbf{{\bar{v}}}}\otimes {\mathbf{{\bar{v}}}}\,+\, \rho \, T  \,{\rm\bf I}{\aftergroup\egroup\originalright}) \, =\, \bm{0}_{{\mathbb{R}}^{D_v}},                   \\
        &\, \\
        & \partial_t E + \operatorname{div}_{\mathbf{x}}{\mathopen{}\mathclose\bgroup\originalleft}( {\mathbf{{\bar{v}}}}{\mathopen{}\mathclose\bgroup\originalleft}( E +\rho \, T{\aftergroup\egroup\originalright})  {\aftergroup\egroup\originalright}) \,=\, 0,
    \end{aligned}
    {\aftergroup\egroup\originalright}.$$ in which $E$ is the second moment of ${f^{\varepsilon}}$, namely the total energy of the gas: $$E = \int_{{\mathbb{R}}^{D_v}} \vert {\mathbf{{\bar{v}}}}\vert^2 f({\mathbf{v}}) \,d{\mathbf{v}}.$$

In the following, we will present the two main collisional kinetic equations that we will consider in the remainder of this paper: the Boltzmann equation (Section \[subsec:boltzmann\_equation\]) and the BGK equation (Section \[subsec:bgk\_equation\]).

Boltzmann equation {#subsec:boltzmann_equation}
------------------

The Boltzmann equation constitutes the cornerstone of the kinetic theory of rarefied gases [@Villani:2002handbook; @CIP]. In a dimensionless, scalar setting, it describes the evolution of the one-particle mass distribution function ${f^{\varepsilon}}({\mathbf{x}},{\mathbf{v}},t) \in {\mathbb{R}}^{+}$, solution to the model equation , in which we still need to specify the collision operator ${{\mathcal{Q}}{\mathopen{}\mathclose\bgroup\originalleft}({f^{\varepsilon}}{\aftergroup\egroup\originalright})({\mathbf{v}})}$. The Boltzmann collision operator models binary elastic collisions between particles having pre-collisional velocities $({\mathbf{v}}',{{\mathbf{v}}_{*}}')$ and post-collisional velocities $({\mathbf{v}},{{\mathbf{v}}_{*}})$. In a two-dimensional velocity space, the pre- and post-collisional velocities are linked through the following parametrization: $${\mathbf{v}}' = \dfrac{{\mathbf{v}}+ {{\mathbf{v}}_{*}}}{2} + \dfrac{{{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}- {{\mathbf{v}}_{*}}{\aftergroup\egroup\originalright}\vert}}{2}{\boldsymbol{\sigma}}, \qquad\quad
    {{\mathbf{v}}_{*}}' = \dfrac{{\mathbf{v}}+ {{\mathbf{v}}_{*}}}{2} - \dfrac{{{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}- {{\mathbf{v}}_{*}}{\aftergroup\egroup\originalright}\vert}}{2}{\boldsymbol{\sigma}},$$ where ${\boldsymbol{\sigma}}$ is the unit vector on the unit circle ${\mathbb{S}}^1 = \{{\boldsymbol{\sigma}}\in {\mathbb{R}}^2: {{\mathopen{}\mathclose\bgroup\originalleft}\vert{\boldsymbol{\sigma}}{\aftergroup\egroup\originalright}\vert} = 1\}$ directed along the pre-collisional relative velocity ${{\mathbf{v}}_r}' = {\mathbf{v}}' - {{\mathbf{v}}_{*}}'$: $${\boldsymbol{\sigma}}= \frac{{{\mathbf{v}}_r}'}{{{\mathopen{}\mathclose\bgroup\originalleft}\vert{{\mathbf{v}}_r}'{\aftergroup\egroup\originalright}\vert}} = \frac{{\mathbf{v}}' - {{\mathbf{v}}_{*}}'}{{{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}' - {{\mathbf{v}}_{*}}'{\aftergroup\egroup\originalright}\vert}}.$$ The Boltzmann collision operator then reads: $$\begin{aligned}
    {{\mathcal{Q}}{\mathopen{}\mathclose\bgroup\originalleft}({f^{\varepsilon}}{\aftergroup\egroup\originalright})({\mathbf{v}})} &= \int_{{\mathbb{R}}^2}\int_{{\mathbb{S}}^1} B({{\mathopen{}\mathclose\bgroup\originalleft}\vert{{\mathbf{v}}_r}{\aftergroup\egroup\originalright}\vert},{\boldsymbol{\sigma}})(f'f_*' - ff_*) d{\boldsymbol{\sigma}}d{{\mathbf{v}}_{*}}\notag \\
    &= \int_{{\mathbb{R}}^2}\int_{0}^{2\pi} B({{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}-{\mathbf{v}}_*{\aftergroup\egroup\originalright}\vert},{\theta_\sigma})(f'f_*' - ff_*) d{\theta_\sigma}d{{\mathbf{v}}_{*}}, \label{eq:collision_operator}\end{aligned}$$ where ${\theta_\sigma}$ is the angle between ${{\mathbf{v}}_r}'$ and ${\boldsymbol{\sigma}}$ and we used the shorthand notations $f = {f^{\varepsilon}}({\mathbf{v}})$, $f_* = {f^{\varepsilon}}({\mathbf{v}}_*)$, $f^{'} = {f^{\varepsilon}}({\mathbf{v}}')$, and $f_*^{'} = {f^{\varepsilon}}({\mathbf{v}}_* ^{'})$. Furthermore, the non-negative function $B({{\mathopen{}\mathclose\bgroup\originalleft}\vert{{\mathbf{v}}_r}{\aftergroup\egroup\originalright}\vert},{\boldsymbol{\sigma}}) \equiv B({{\mathopen{}\mathclose\bgroup\originalleft}\vert{{\mathbf{v}}_r}{\aftergroup\egroup\originalright}\vert},{\theta_\sigma})$ is the collision kernel, which, by physical arguments of invariance, only depends on the relative speed ${{\mathopen{}\mathclose\bgroup\originalleft}\vert{{\mathbf{v}}_r}{\aftergroup\egroup\originalright}\vert} = {{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}-{\mathbf{v}}_*{\aftergroup\egroup\originalright}\vert}$ and $\cos({\theta_\sigma}) = {{\mathbf{v}}_r}'/{{\mathopen{}\mathclose\bgroup\originalleft}\vert{{\mathbf{v}}_r}{\aftergroup\egroup\originalright}\vert} \cdot {\boldsymbol{\sigma}}$. The collision kernel $B$ contains all relevant microscopic information such as the kind of particles and type of interactions. For instance, when particles interact via an inverse power law potential ${\Phi}(r) = r^{-k+1}$ $(k > 2)$, with $r$ the inter-particle distance, $B$ factors as: $$\label{eq:collision_kernel_factors_sigma}
    B({{\mathopen{}\mathclose\bgroup\originalleft}\vert{{\mathbf{v}}_r}{\aftergroup\egroup\originalright}\vert},{\theta_\sigma}) = {{\mathopen{}\mathclose\bgroup\originalleft}\vert{{\mathbf{v}}_r}{\aftergroup\egroup\originalright}\vert}^\gamma b_\gamma({\theta_\sigma}), \qquad\quad
    \gamma = \frac{k-3}{k-1}.$$ Notably, in the special case $k = 3$, the collision kernel in is independent of the relative speed ${{\mathopen{}\mathclose\bgroup\originalleft}\vert{{\mathbf{v}}_r}{\aftergroup\egroup\originalright}\vert}$ and the resulting particles are known as Maxwellian particles. If, in addition, $b_0({\theta_\sigma}) = b_0$ is assumed constant, the particles are referred to as pseudo-Maxwellian particles. In general (except for hard sphere and pseudo-Maxwellian particles), the angular collision kernel $b_\gamma$ in is expressed implicitly and contains a singularity for grazing collisions $({\theta_\sigma}\to 0)$, and its mathematical analysis, as well as its numerical simulations, can be very difficult [@alexandre2000entropy]. For that reason, the angular collision kernel is usually replaced by an integrable function by cutting off such grazing collision angles (Grad’s cut-off assumption) [@Cercignani1988].

It is instructive to split the collision operator into a gain and loss operator as: $$\label{eq:collision_operator_gain_loss}
    {{\mathcal{Q}}{\mathopen{}\mathclose\bgroup\originalleft}({f^{\varepsilon}}{\aftergroup\egroup\originalright})({\mathbf{v}})} = {{{\mathcal{Q}}^+}{\mathopen{}\mathclose\bgroup\originalleft}({f^{\varepsilon}}{\aftergroup\egroup\originalright})({\mathbf{v}})} - {\nu}({f^{\varepsilon}}){f^{\varepsilon}}({\mathbf{v}}),$$ where the gain operator is given by: $$\label{eq:collision_operator_gain}
    {{{\mathcal{Q}}^+}{\mathopen{}\mathclose\bgroup\originalleft}({f^{\varepsilon}}{\aftergroup\egroup\originalright})({\mathbf{v}})} = \int_{{\mathbb{R}}^2}\int_{0}^{2\pi} B({{\mathopen{}\mathclose\bgroup\originalleft}\vert{{\mathbf{v}}_r}{\aftergroup\egroup\originalright}\vert},{\theta_\sigma})f'f_*' d{\theta_\sigma}d{{\mathbf{v}}_{*}}.$$ The loss operator ${{\mathcal{Q}}^-}({f^{\varepsilon}}) = {\nu}({f^{\varepsilon}}){f^{\varepsilon}}$ contains the collision frequency ${\nu}({f^{\varepsilon}}) \in {\mathbb{R}}^+$, which is defined by: $$\label{eq:collision_frequency_Boltzmann}
    {\nu}({f^{\varepsilon}}) = \int_{{\mathbb{R}}^2}\int_{0}^{2\pi} B({{\mathopen{}\mathclose\bgroup\originalleft}\vert{{\mathbf{v}}_r}{\aftergroup\egroup\originalright}\vert},{\theta_\sigma})f_* d{\theta_\sigma}d{{\mathbf{v}}_{*}}.$$ (Evidently, equation is only valid if both integrals - are convergent, which is certainly true for a cut-off collision kernel.) In general, the collision frequency depends on the dimension of velocity space and the type of microscopic collisions. In particular, when considering pseudo-Maxwellian particles for ${D_v = 2}$, implying that $\gamma = 0$ and $b_0({\theta_\sigma}) = b_0$ constant, the collision kernel becomes much simpler and is given by ${B({{\mathopen{}\mathclose\bgroup\originalleft}\vert{{\mathbf{v}}_r}{\aftergroup\egroup\originalright}\vert},{\theta_\sigma}) = b_0}$, that is, independent of ${{\mathopen{}\mathclose\bgroup\originalleft}\vert{{\mathbf{v}}_r}{\aftergroup\egroup\originalright}\vert}$ and ${\theta_\sigma}$. In that case, the collision frequency can be explicitly computed as: $${\nu}({f^{\varepsilon}}) = b_0 \int_{0}^{2\pi} d{\theta_\sigma}\int_{{\mathbb{R}}^2} f_* d{{\mathbf{v}}_{*}}= 2\pi b_0\rho,$$ and the collision operator in reads: $$\label{eq:collision_operator_maxwellian_particles}
    {{\mathcal{Q}}{\mathopen{}\mathclose\bgroup\originalleft}({f^{\varepsilon}}{\aftergroup\egroup\originalright})({\mathbf{v}})} = {{{\mathcal{Q}}^+}{\mathopen{}\mathclose\bgroup\originalleft}({f^{\varepsilon}}{\aftergroup\egroup\originalright})({\mathbf{v}})} - 2\pi b_0\rho {f^{\varepsilon}}.$$ We note that, for inverse power law potentials with $\gamma \neq 0$ (including hard sphere particles), the collision frequency in general depends on the density, temperature and collision kernel, see [@Struchtrup2005]. Finally, it is easy (see [@CIP]) to check that the Boltzmann collision operator satisfies the hypotheses \[hypConservations\]-\[hypEntropy\]-\[hypEquilib\].

BGK equation {#subsec:bgk_equation}
------------

Due to the high-dimensional and complicated structure of the Boltzmann collision operator, the Boltzmann collision kernel  is often replaced by simpler collision models that capture most of its essential features. The most well-known approximation is the BGK model [@Bhatnagar1954], which models collisions as a relaxation towards thermodynamic equilibrium. (Because of the moments-dependency of the equilibrium, this is still a nonlinear operator.) It is the most well regarded simplified model of the Boltzmann equation, and is almost universally used in the physics and numerics communities (see [@DimarcoPareschi2014] and the references therein for details). It is given by: $$\label{eq:bgk_equation}
    \partial_t {f^{\varepsilon}}+ {\mathbf{v}}\cdot \nabla_{{\mathbf{x}}} {f^{\varepsilon}}= \frac{\nu}{{\varepsilon}}({{{\mathcal{M}}_{{\mathbf{v}}}{\mathopen{}\mathclose\bgroup\originalleft}({f^{\varepsilon}}{\aftergroup\egroup\originalright})}}- {f^{\varepsilon}}),$$ where ${{{\mathcal{M}}_{{\mathbf{v}}}{\mathopen{}\mathclose\bgroup\originalleft}({f^{\varepsilon}}{\aftergroup\egroup\originalright})}}$ denotes the local Maxwellian distribution, which, for a $D_v$-dimensional velocity space, is given by: $$\label{eq:maxwellian} 
    {{{\mathcal{M}}_{{\mathbf{v}}}{\mathopen{}\mathclose\bgroup\originalleft}({f^{\varepsilon}}{\aftergroup\egroup\originalright})}}= \frac{\rho}{(2\pi T)^{D_v/2}} \exp{{\mathopen{}\mathclose\bgroup\originalleft}(-\frac{|{\mathbf{v}}-{\mathbf{{\bar{v}}}}|^2}{2T}{\aftergroup\egroup\originalright})} := {\mathcal{M}}_{\mathbf{v}}^{\rho,{\mathbf{{\bar{v}}}},T}.$$ Furthermore, in the BGK equation , the collision frequency ${\nu}\in {\mathbb{R}}^+$ is in general derived from the Boltzmann collision operator and its expression is given in , while in the linearized setting, ${\nu}$ is independent of ${f^{\varepsilon}}$ and is formulated in , see Section \[subsec:lin\_boltzmann\_equation\]. Notably, when setting ${\nu}= \rho$ for $D_v = 2$, the BGK model matches the loss term of the Boltzmann collision operator for pseudo-Maxwellian particles close to equilibrium, see equation . The BGK collision operator satisfies by construction the hypotheses \[hypConservations\]-\[hypEntropy\]-\[hypEquilib\].

Numerical method {#sec:time_integrator}
================

Now that we have introduced the model problems, we turn to the description of the numerical method that will be the focus of this paper. Equation  needs to be discretized in space, velocity and time.

We discretize equation in space using finite differences on a uniform, constant in time, periodic mesh with spacing ${\Delta x}$, consisting of $I$ mesh points $x_i=i{\Delta x}$, ${1 \le i \le I}$, with $I{\Delta x}=1$. In the numerical experiments of section \[sec:results\], we use the classical WENO scheme [@shu:1998], which we briefly recall in Appendix \[sub:weno\]. Next, we discretize velocity space by choosing $J$ discrete components denoted by ${\mathbf{v}}_j$. For the Boltzmann equation, that is, equation with collision operator , we use the fast spectral discretization of the Boltzmann operator, taken from [@MoPa:2006]. This method is recalled in Appendix \[sub:FastSpectral\].

The semidiscrete numerical solution on this mesh is denoted by $f_{i,j}(t)$, where we have dropped the superscript ${\varepsilon}$ on discretized quantities. We then obtain a semidiscrete system of ODEs of the form: $$\label{eq:semidiscrete} 
    \dot{{\mathbf{f}}} = {\mathrm{D}_t}({\mathbf{f}}),  \qquad
    {\mathrm{D}_t}({\mathbf{f}}) = -{\mathrm{D}_{{\boldsymbol{x}},{\boldsymbol{v}}}}({\mathbf{f}}) + \frac{1}{{\varepsilon}}{\mathcal{Q}}{\mathopen{}\mathclose\bgroup\originalleft}({\mathbf{f}}{\aftergroup\egroup\originalright}),$$ where ${\mathrm{D}_{{\boldsymbol{x}},{\boldsymbol{v}}}}(\cdot)$ represents the finite difference discretization of the convective derivative ${\mathbf{v}}\cdot \nabla_{\mathbf{x}}$, and ${\mathbf{f}}$ is a vector of size $I \cdot J$.

In the remainder of this section, we describe the time discretization of the semi-discretized system , which is the novel element in the full discretization of equation . We start in Section \[sub:projective\_integration\] with the projective integration method, which aims at efficiently simulating systems with *exactly two* time scales (one fast and one slow). In Section \[sub:telescopic\_integration\], we present the generalized telescopic projective integration method, which can deal with multiple fast time scales.

Projective integration {#sub:projective_integration}
----------------------

Projective integration [@Gear2003projective] is a method that is tailored to problems with exactly two distinct time scales. As such, in the context of kinetic equations, it matches nicely with the spectral properties of a linear BGK equation, as was shown in [@Lafitte2012]. Projective integration combines a few small time steps with a naive (*inner*) timestepping method (here, a direct forward Euler discretization) with a much larger (*projective, outer*) time step. The idea is sketched in figure \[fig:proj\_int\].

[![\[fig:proj\_int\] Sketch of projective integration. At each time, an explicit method is applied over a number of small time steps (black dots) so as to stably integrate the fast modes. As soon as these modes are sufficiently damped, the solution is extrapolated using a much larger time step (dashed lines). ](image/PI_sketch "fig:")]{}

#### Inner integrators.

At the innermost leve, we introduce a uniform time mesh with time step ${\delta t}$ and discrete time instants $t^k=k{\delta t}$. At this leve, we choose the (explicit) forward Euler method with time step ${\delta t}$, for which we will, later on, use the shorthand notation: $$\label{eq:fe_scheme} 
    {\mathbf{f}}^{k+1} = {S_{{\delta t}}}({\mathbf{f}}^{k}) = {\mathbf{f}}^k + {\delta t}{\mathrm{D}_t}({\mathbf{f}}^k), \qquad k = 0, 1, \ldots.$$ The purpose of the inner integrator is to capture the fastest components in the numerical solution of system and to sufficiently damp these out. We only require the innermost integrator to be stable for these components. The size of the inner time step ${\delta t}$ and the required number of inner steps $K$ will depend on the spectral properties of the semidiscretization . This will be studied in Section \[sec:linearOpSpecProp\].

#### Outer integrators.

In system , the small parameter ${\varepsilon}$ leads to a classical time step restriction of the form ${\delta t}= O({\varepsilon})$ for the inner integrator. However, as ${\varepsilon}$ goes to $0$, we obtain the limiting system , for which a standard finite volume/forward Euler method only needs to satisfy a CFL stability restriction of the form ${\Delta t}\le C{\Delta x}$, with $C$ a constant that depends on the specific choice of the scheme.

In [@Lafitte2012], it was proposed to use a projective integration method to accelerate such a brute-force integration; the idea, originating from [@Gear2003projective], is the following. Starting from a computed numerical solution ${\mathbf{f}}^n$ at time $t^n=n{\Delta t}$, one first takes $K+1$ *inner* steps of size ${\delta t}$ using , denoted as ${\mathbf{f}}^{n,k+1}$, in which the superscripts $(n,k)$ denote the numerical solution at time ${t^{n,k}=n{\Delta t}+k{\delta t}}$. The aim is to obtain a discrete derivative to be used in the *outer* step to compute ${\mathbf{f}}^{n+1} = {\mathbf{f}}^{n+1,0}$ via extrapolation in time: $$\begin{aligned}
    {\mathbf{f}}^{n+1} & = {\mathbf{f}}^{n,K+1} + ({\Delta t}- (K + 1){\delta t})\frac{{\mathbf{f}}^{n,K+1} - {\mathbf{f}}^{n,K}}{{\delta t}}, \\
     & = {\mathbf{f}}^{n,K+1} + M{\delta t}\frac{{\mathbf{f}}^{n,K+1} - {\mathbf{f}}^{n,K}}{{\delta t}},\end{aligned}$$ where $M = {\Delta t}/{\delta t}-(K+1)$. Also the size of the (macroscopic) extrapolation step ${\Delta t}$ will result from the spectral analysis of the semidiscretization  in section \[sec:linearOpSpecProp\].

Higher-order projective Runge-Kutta (PRK) methods have been constructed [@LafitteLejonSamaey2015; @LafitteMelisSamaey2017] by replacing each time derivative evaluation $\mathbf{k}_s$ in a classical Runge-Kutta method by $K+1$ steps of an inner integrator as follows: $$\begin{aligned}
    s = 1 :\;\; & 
    \begin{dcases} 
        {\mathbf{f}}^{n,k+1} &= {\mathbf{f}}^{n,k} + {\delta t}{\mathrm{D}_t}({\mathbf{f}}^{n,k}), \qquad 0 \le k \le K \\ 
        \mathbf{k}_1 &= \dfrac{{\mathbf{f}}^{n,K+1} - {\mathbf{f}}^{n,K}}{{\delta t}}
    \end{dcases} \\
    2 \le s \le S :\;\; & 
    \begin{dcases} 
        {\mathbf{f}}^{n+c_s,0}_s &= {\mathbf{f}}^{n,K+1} + (c_s{\Delta t}-(K+1){\delta t}) \sum_{l=1}^{s-1}\dfrac{a_{s,l}}{c_s} \mathbf{k}_l, \\
        {\mathbf{f}}^{n+c_s,k+1}_s &= {\mathbf{f}}^{n+c_s,k}_s + {\delta t}{\mathrm{D}_t}({\mathbf{f}}^{n+c_s,k}_s), \qquad 0 \le k \le K \\
        \mathbf{k}_s &= \dfrac{{\mathbf{f}}^{n+c_s,K+1}_s - {\mathbf{f}}^{n+c_s,K}_s}{{\delta t}}
    \end{dcases} \\
    & {\mathbf{f}}^{n+1} = {\mathbf{f}}^{n,K+1} + ({\Delta t}-(K+1){\delta t})\sum_{s=1}^{S}b_s \mathbf{k}_s.\end{aligned}$$ To ensure consistency, the Runge-Kutta matrix $\mathbf{a}=(a_{s,i})_{s,i=1}^S$, weights ${\mathbf{b}=(b_s)_{s=1}^S}$, and nodes $\mathbf{c}=(c_s)_{s=1}^S$ satisfy the conditions $0\le b_s \le 1$ and $0 \le c_s \le 1,$ as well as: $$\sum_{s=1}^Sb_s=1, \qquad \sum_{i=1}^{S-1} a_{s,i} =c_s, \quad 1 \le s \le S.$$

Telescopic projective integration {#sub:telescopic_integration}
---------------------------------

In general, the stiff semidiscrete system , contains more than two distinct time scales. In this section, we therefore describe an extension of projective integration, called *telescopic projective integration* (TPI) and introduced in [@Gear2003telescopic], that can handle multiple time scales. This method has been studied in the context of linear BGK equations with multiple relaxation times in [@MelisSamaey2017].

Telescopic projective integration employs a number of projective integrator levels, which, starting from a base (*innermost*) integrator, are wrapped around the previous level integrator [@Gear2003telescopic]. In this way, a hierarchy of projective integrators is formed in which each level (except the innermost and outermost one) fulfils both an inner and outer integrator role. This generalizes the idea of projective integration, which contains only one projective level wrapped around an inner integrator. The idea of a level-3 TPI method with $K=2$ on each projective level is sketched in figure \[fig:tpi\_sketch\]. The different level integrators in a TPI method can in principle be selected independently from each other, but in general one selects a first order explicit scheme (the forward Euler scheme) for all but the outermost integrator level, whose order is chosen to meet the accuracy requirements dictated by the problem.

[![\[fig:tpi\_sketch\] A level-3 TPI method drawn for three outermost time steps $h_3$ (bottom row) with $K=2$ on all projective levels. The dots correspond to different time points at which the numerical solution is calculated. The time step and projective step size of each level $\ell=0,\ldots,2$ are denoted by $h_\ell$ and $M_\ell$, respectively. ](image/TPI_level3_sketch "fig:")]{}

#### Innermost integrator

We intend to integrate the semidiscrete system of equations using a uniform time mesh with time step $h_0$ and discrete time instants $t^k=kh_0$. The innermost integrator of the TPI method is chosen to be the forward Euler (FE) method, $${\mathbf{f}}^{k+1} = {\mathbf{f}}^k + h_0{\mathrm{D}_t}({\mathbf{f}}^k).$$ In the sequel, we use the following shorthand notation: $${\mathbf{f}}^{k+1} = S_{0}({\mathbf{f}}^{k}) \qquad (k = 0, 1, \ldots),$$ in which $S_0$ denotes the time stepper with corresponding time step $h_0$. Also in the telescopic projective integration method, the purpose of the innermost integrator is only to capture the fastest components in the numerical solution of system  and to sufficiently damp these out. As a consequence, it is ill-advised to use higher-order methods for the innermost integrator, see [@MelisSamaey2017] for a more detailed discussion.

#### Projective (outer) levels

The telescopic projective integration method employs in general $L$ nested levels of projective integration that are constructed around the innermost integrator. In [@Gear2003telescopic], the method has been introduced in a recursive way. Here, following [@MelisSamaey2017], we describe the method in an alternative way, to make the presentation more similar to that of classical projective integration.

To keep track of the time instant at which the numerical solution is computed throughout the TPI method and at the same time desiring a compact notation, in what follows, we employ superscript triplets of the form $(\ell,n,k_\ell)$ where $\ell$ denotes the integrator level ranging from $0$ (innermost) to $L-1$, $n$ represents the index of the current outermost integrator time $t^n=nh_L$, and $k_\ell$ corresponds to the iteration index of the integrator on level $\ell$. The numerical time on each level $\ell = 0, \ldots, L-1$ is then defined as (see also figure \[fig:tpi\_sketch\]): $$\label{eq:tpi_time}
    t^{\ell,n,k_\ell} = nh_L + \sum_{\ell'=\ell}^{L-1} k_{\ell'}h_{\ell'}.$$ Notice that, for a certain level $\ell$, this time requires the iteration indices $k_{\ell'}$ of all its outer integrators. Therefore, it incorporates a memory that keeps up with the current time instants at which the outer integrators of a given level $\ell$ integrator have arrived at and is necessary to take into account to correctly reflect the numerical time of the solution on each level $\ell$.

Starting from a computed numerical solution ${\mathbf{f}}^n$ at time $t^n=nh_L$, one first takes $K_0+1$ steps of size $h_0$ with the innermost integrator: $$\label{eq:semidiscrete_innermost_scheme}
    {\mathbf{f}}^{0,n,k_0+1} =S_{0}({\mathbf{f}}^{0,n,k_0}) \qquad (0 \le k_0 \le K_0),$$ in which ${\mathbf{f}}^{0,n,k_0}$ corresponds to the numerical solution at time $t^{0,n,k_0}$ calculated by the innermost integrator. Since all outer integrator iteration indices $k_{\ell'}$, $\ell' = 1,\ldots,L-1$ are initially zero in , we have ${t^{0,n,k_0} = nh_L+k_0h_0}$. The repeated action of the innermost integrator is depicted by small black arrows in the upper row of figure \[fig:tpi\_sketch\], for which we chose $K_0=2$.

In the telescopic projective integration framework, the scheme is set up from the lowest level up to the highest level. The aim is to obtain a discrete derivative to be used on each level to eventually compute ${\mathbf{f}}^{n+1} = {\mathbf{f}}^{0,n+1,0}$ via extrapolation in time. Using the innermost integrator iterations , we perform the extrapolation by a projective integrator on level 1, written as: $$\label{eq:tpfe_first_level} 
    {\mathbf{f}}^{1,n,1} = {\mathbf{f}}^{0,n,K_{0}+1} + {\mathopen{}\mathclose\bgroup\originalleft}(M_{0}h_{0}{\aftergroup\egroup\originalright})\frac{{\mathbf{f}}^{0,n,K_{0}+1} - {\mathbf{f}}^{0,n,K_{0}}}{h_{0}},$$ which corresponds to the projective forward Euler (PFE) method. In , ${\mathbf{f}}^{1,n,1}$ represents the numerical solution at time $t^{1,n,1}$ calculated by one iteration of the first level projective integrator. Since $k_1=1$ and all its outer integrator iteration indices $k_{\ell'}$, ${\ell'=2,\ldots,L-1}$ are still zero in , we have ${t^{1,n,1} = nh_L + h_1}$. One step of the first level integrator is visualized by a large green arrow in the upper row of figure \[fig:tpi\_sketch\]. By repeating this idea, we construct a hierarchy of projective integrators on levels $\ell=1,\ldots,L-1$, given by: $$\label{eq:tpfe_projective_levels} 
    {\mathbf{f}}^{\ell,n,k_\ell+1} = {\mathbf{f}}^{\ell-1,n,K_{\ell-1}+1} + {\mathopen{}\mathclose\bgroup\originalleft}(M_{\ell-1}h_{\ell-1}{\aftergroup\egroup\originalright})\frac{{\mathbf{f}}^{\ell-1,n,K_{\ell-1}+1} - {\mathbf{f}}^{\ell-1,n,K_{\ell-1}}}{h_{\ell-1}},$$ in which, on each level $\ell$, we iterate over $k_{\ell} = 0, \ldots, K_{\ell}$. In , ${\mathbf{f}}^{\ell,n,k_\ell}$ denotes the numerical solution at time $t^{\ell,n,k_\ell}$ calculated by the projective integrator on level $\ell$. According to , this time depends on the values $k_{\ell'}$, $\ell'=\ell+1,\ldots,L-1$ of all of its outer integrators. In figure \[fig:tpi\_sketch\], these projective integrator steps are shown by long arrows for each level $\ell=1,\ldots,3$. Ultimately, the outermost integrator on level $L$ computes ${\mathbf{f}}^{n+1}$ as: $$\label{eq:tpfe_outermost}
    {\mathbf{f}}^{n+1} = {\mathbf{f}}^{L-1,n,K_{L-1}+1} + (M_{L-1}h_{L-1})\frac{{\mathbf{f}}^{L-1,n,K_{L-1}+1} - {\mathbf{f}}^{L-1,n,K_{L-1}}}{h_{L-1}}.$$ Since the outermost integrator also constitutes a PFE scheme, the telescopic method resulting from the hierarchy of projective levels - is called telescopic projective forward Euler (TPFE).

It is straightforward to implement higher-order extensions of the outermost integrator, as is done in [@MelisSamaey2017]. We mention the projective Runge-Kutta methods of order 2 and 4, leading to TPRK2 and TPRK4 method in the telescopic case. In general, the outermost integrator in a TPRK method replaces each time derivative evaluation $\mathbf{k}_s$ in a classical Runge-Kutta method by $K_{L-1}+1$ steps of its inner integrator on level $L-1$. Using with $\ell=L-1$, the first stage in a TPRK method calculates the time derivative $\mathbf{k}_1$ as: $$\label{eq:tprk_stage_1}
    \mathbf{k}_1 = \dfrac{{\mathbf{f}}^{L-1,n,K_{L-1}+1} - {\mathbf{f}}^{L-1,n,K_{L-1}}}{h_{L-1}}.$$ Computing $\mathbf{k}_s$ on any other stage $s \ge 2$ requires evaluating time derivatives at the intermediate times ${t^{n+c_s} = (n+c_s)h_L}$. Similarly to , these are computed as: $$\label{eq:tprk_stage_s}
    \mathbf{k}_s = \dfrac{{\mathbf{f}}^{L-1,n+c_s,K_{L-1}+1} - {\mathbf{f}}^{L-1,n+c_s,K_{L-1}}}{h_{L-1}}.$$ However, since the numerical solution at time $t^{n+c_s}$ in equation is not available, we use the integrator on level $L-1$ to approximate it as follows: $$\begin{dcases}
        {\mathbf{f}}^{L-1,n+c_s,0} = {\mathbf{f}}^{L-1,n,K_{L-1}+1} + (c_sh_L-(K_{L-1}+1)h_{L-1}) \sum_{i=1}^{s-1}\dfrac{a_{s,i}}{c_s} \mathbf{k}_i \\
        {\mathbf{f}}^{L-1,n+c_s,k_{L-1}+1} = {\mathbf{f}}^{L-2,n+c_s,K_{L-2}+1} \\
        \hspace{3.1cm} + \,{\mathopen{}\mathclose\bgroup\originalleft}(M_{L-2}h_{L-2}{\aftergroup\egroup\originalright})\frac{{\mathbf{f}}^{L-2,n+c_s,K_{L-2}+1} - {\mathbf{f}}^{L-2,n+c_s,K_{L-2}}}{h_{L-2}},
    \end{dcases}$$ in which the second equation iterates over $0 \le k_{L-1} \le K_{L-1}$. Ultimately, the outermost integrator of a TPRK method is written as: $${\mathbf{f}}^{n+1} = {\mathbf{f}}^{L-1,n,K_{L-1}+1} + (M_{L-1}h_{L-1})\sum_{s=1}^{S}b_s \mathbf{k}_s.$$ The consistency conditions on the Runge-Kutta matrix $\mathbf{a}=(a_{s,i})_{s,i=1}^S$, weights $\mathbf{b}=(b_s)_{s=1}^S$, and nodes ${\mathbf{c}=(c_s)_{s=1}^S}$ are still valid in this setting [@MelisSamaey2017]. In the numerical experiments in Section \[sec:results\], we will use the projective Runge-Kutta method of order 4 as outermost integrator.

On linearized operators and spectral properties {#sec:linearOpSpecProp}
===============================================

Telescopic projective integration methods are very versatile, but require choosing a relatively large number of method parameters: the size of the time steps $h_\ell$ at each level, the number $K_\ell$ of inner steps at each level, as well as the number $L$ of telescopic levels. As these choices are dictated mainly by stability requirements, they crucially depend on the spectrum of the collision operator. The analysis of this spectrum for the problems of Section \[sec:models\] is the focus of this Section. In [@Lafitte2012; @LafitteLejonSamaey2015], the spectrum of the collision operator was analysed for linear kinetic equations in the diffusive and hyperbolic scalings. An extension to kinetic relaxations of a nonlinear hyperbolic conservation law was presented in [@LafitteMelisSamaey2017]. In all of these settings, the spectrum of the collision operator turned out to consist of exactly two well-separated time scales, and projective integration was therefore sufficient. A first study of telescopic projective integration for linear kinetic equations with multiple relaxation times was presented in [@MelisSamaey2017].

In this section, we devise a general framework of linear operators in which linearizations of both the BGK and Boltzmann equation can be studied. This framework will allow determining suitable method parameters for the (telescopic) projective integration of the Boltzmann and nonlinear BGK equations. In addition, it allows embedding the linear kinetic equations that were studied in [@LafitteMelisSamaey2017; @MelisSamaey2017]. For the reader’s convenience, we restrict the exposition to the case when $D_v = 2$. However, all the results of this section can be extended straigthforwardly to the $D_v=3$ case, at the cost of heavier notations. In Section \[subsec:lin\_bgk\], we build this framework for BGK equations and draw conclusions on their spectral properties. Afterwards, in Section \[subsec:lin\_boltzmann\_equation\], we extend this framework to include the Boltzmann equation. We discuss the selection of suitable method parameters for (telescopic) projective integration in Section \[sec:method\_param\].

Linearized BGK models and their spectra {#subsec:lin_bgk}
---------------------------------------

We first recall the linearization of the BGK equation , as it has been described in [@Cercignani1988]. We then show in section \[sec:linear\_kinetic\] how the simpler linear kinetic equations that were analyzed in [@LafitteMelisSamaey2017; @MelisSamaey2017] fit in this framework. Finally, in section \[sec:bgk\_spectrum\], we discuss how the analysis of the spectrum of linear kinetic equations in [@LafitteMelisSamaey2017; @MelisSamaey2017] generalizes to the linearization of the full BGK equation.

### Linearized BGK equation {#sec:linearized_bgk}

In [@Cercignani1988], it is shown that the linearized BGK operator can be formulated as: $$\label{eq:linearized_bgk_maxwellian}
    {\mathcal{M}}_1({f^{\varepsilon}})({\mathbf{x}},{\mathbf{v}},t) = \sum_{k=0}^{D_v+1} \Psi_k({\mathbf{v}}) \langle\Psi_k, {f^{\varepsilon}}\rangle({\mathbf{x}},t),$$ in which the scalar product is defined by: $$\label{eq:scalar_product_Hilbert}
    \langle g,h \rangle= \int_{{\mathbb{R}}^{D_v}} g({\mathbf{v}})\overline{h({\mathbf{v}})} \frac{1}{2\pi}\exp{\mathopen{}\mathclose\bgroup\originalleft}(\frac{-{{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}{\aftergroup\egroup\originalright}\vert}^2}{2}{\aftergroup\egroup\originalright}) d{\mathbf{v}}.$$ Furthermore, the basis functions $\Psi_k({\mathbf{v}})$ in represent an orthogonal basis for the space spanned by $$V:= {\mathopen{}\mathclose\bgroup\originalleft}\{1, v^{x_1}, \ldots, v^{x_{D_v}}, \frac{{{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}{\aftergroup\egroup\originalright}\vert}^2}{2}{\aftergroup\egroup\originalright}\},$$ in which the superscript $x_d$ on $v$ indicates that $v^{x_d}$ is the component of ${\mathbf{v}}$ in the $d$-th velocity dimension. The space $V$ corresponds to the set of elementary collision invariants $(\psi_k({\mathbf{v}}))_{k=0}^{{D_v}+1}$. We construct an orthonormal basis, i.e., we seek basis functions $Psi_k$ such that: $$\label{eq:normalization_Psi}
    (\Psi_k, \Psi_j) = \delta_{kj}, \qquad k,j \in \{0,\ldots,{D_v}+1\},$$ in which $\delta_{kj}$ denotes the Kronecker delta. A straightforward application of the Gram-Schmidt process then allows to compute the desired set of orthonormal functions $\Psi_k({\mathbf{v}})$ satisfying as: $$\label{eq:Psi_normalized}
    \big(\Psi_0({\mathbf{v}}), \ldots, \Psi_{D_v+1}({\mathbf{v}})\big) = {\mathopen{}\mathclose\bgroup\originalleft}(1, v^{x_1},\ldots, v^{x_{D_v}}, \frac{{{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}{\aftergroup\egroup\originalright}\vert}^2 - 2}{2}{\aftergroup\egroup\originalright}).$$ Ultimately, using the linearized Maxwellian , the linearized version of the full BGK equation reads: $$\label{eq:linearized_bgk_equation_framework}
    \partial_t {f^{\varepsilon}}+ {\mathbf{v}}\cdot \nabla_{{\mathbf{x}}} {f^{\varepsilon}}= -\frac{{\nu}}{{\varepsilon}}({\mathcal{I}}- \Pi_\text{BGK}){f^{\varepsilon}},$$ where $\Pi_\text{BGK}$ is the following rank-${D_v}+2$ projection operator: $$\label{eq:projection_operator_bgk}
    \Pi_\text{BGK} {f^{\varepsilon}}= \sum_{k=0}^{{D_v}+1} \Psi_k({\mathbf{v}})\langle\Psi_k, {f^{\varepsilon}}\rangle.$$

### Linear kinetic models\[sec:linear\_kinetic\]

In [@LafitteMelisSamaey2017], the spectrum of a specific class of linear, hyperbolically scaled, kinetic equations was studied. In a scalar, two-dimensional setting (where we shall set $x := x_1$ and $y := x_2$), the following equation was proposed: $$\label{eq:artificial_bgk_equation}
    \partial_t {f^{\varepsilon}}+ {\mathbf{v}}\cdot \nabla_{{\mathbf{x}}} {f^{\varepsilon}}= \frac{1}{{\varepsilon}}({{{\mathcal{M}}_{{\mathbf{v}}}{\mathopen{}\mathclose\bgroup\originalleft}(\rho^{\varepsilon}{\aftergroup\egroup\originalright})}} - {f^{\varepsilon}}),$$ with an artificial Maxwellian distribution given by: $$\label{eq:artificial_maxwellian}
    {{{\mathcal{M}}_{{\mathbf{v}}}{\mathopen{}\mathclose\bgroup\originalleft}(\rho^{\varepsilon}{\aftergroup\egroup\originalright})}} = \rho^{\varepsilon}(1+  v^x + v^y),$$ in which $\rho^{\varepsilon}$ is linked to ${f^{\varepsilon}}$ by averaging over the velocity space: $$\rho^{\varepsilon}= \int_{{\mathbb{R}}^2} {f^{\varepsilon}}({\mathbf{v}}) \frac{1}{2\pi} \exp{\mathopen{}\mathclose\bgroup\originalleft}(\frac{-{{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}{\aftergroup\egroup\originalright}\vert}^2}{2}{\aftergroup\egroup\originalright})d{\mathbf{v}}$$ see also equation . In [@MelisSamaey2017], multiple relaxation time linearized BGK equations of the following form are considered: $$\label{eq:artificial_linearized_bgk_equation}
    \partial_t {f^{\varepsilon}}+ {\mathbf{v}}\cdot \nabla_{{\mathbf{x}}} {f^{\varepsilon}}= \frac{{\nu}}{{\varepsilon}}({{{\mathcal{M}}_{\text{\scriptsize{lin}},{\mathbf{v}}}{\mathopen{}\mathclose\bgroup\originalleft}({f^{\varepsilon}}{\aftergroup\egroup\originalright})}}- {f^{\varepsilon}}),$$ with linearized Maxwellian given by: $$\label{eq:artificial_maxwellian_linearized} 
    {{{\mathcal{M}}_{\text{\scriptsize{lin}},{\mathbf{v}}}{\mathopen{}\mathclose\bgroup\originalleft}({f^{\varepsilon}}{\aftergroup\egroup\originalright})}}= \rho^{\varepsilon}(1 + v^x)(1 + v^y).$$

To fit both of these model problems in the general framework of section \[sec:linearized\_bgk\], we rewrite the artificial Maxwellians in and as: $$\begin{aligned}
    {\mathcal{M}}_{-2}({f^{\varepsilon}})({\mathbf{x}},{\mathbf{v}},t) &= \Psi_{-2}({\mathbf{v}}) \rho^{\varepsilon}({\mathbf{x}},t),& \quad \Psi_{-2}({\mathbf{v}}) &= (1 + v^x + v^y), \label{eq:framework_artificial} \\
    {\mathcal{M}}_{-1}({f^{\varepsilon}})({\mathbf{x}},{\mathbf{v}},t) &= \Psi_{-1}({\mathbf{v}}) \rho^{\varepsilon}({\mathbf{x}},t),& \quad \Psi_{-1}({\mathbf{v}}) &= (1 + v^x)(1 + v^y). \label{eq:framework_linearized_artificial}\end{aligned}$$ With these choices of the Maxwellian distribution, we can summarize the linear kinetic models and as: $$\partial_t {f^{\varepsilon}}+ {\mathbf{v}}\cdot \nabla_{{\mathbf{x}}} {f^{\varepsilon}}= -\frac{{\nu}}{{\varepsilon}}({\mathcal{I}}- \Pi_k){f^{\varepsilon}},$$ where ${\mathcal{I}}$ is the identity operator and $\Pi_k$ is the following rank-1 projection operator: $$\label{eq:projection_operator_artificial}
    \Pi_k {f^{\varepsilon}}= \Psi_k({\mathbf{v}}) \int_{{\mathbb{R}}^2} {f^{\varepsilon}}({\mathbf{v}}) \frac{1}{2\pi}\exp{\mathopen{}\mathclose\bgroup\originalleft}(\frac{-{{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}{\aftergroup\egroup\originalright}\vert}^2}{2}{\aftergroup\egroup\originalright}) d{\mathbf{v}},$$ with either $k = -2$ or $k = -1$ .

### Linearized BGK spectrum\[sec:bgk\_spectrum\]

The above considerations indicate that the structure of the linearized Maxwellian and the linearized BGK projection operator are almost identical to those in - and , respectively. Indeed, the linear kinetic projection operators $\Pi_{-2}$ and $\Pi_{-1}$ in match the first three terms of $\Pi_\text{BGK}$ in and only differ in the last term. This can be seen by using the orthonormal set of basis functions and the scalar product , and subsequently rewriting $\Pi_{-2}$ and $\Pi_{-1}$ as: $$\label{eq:projection_operator_artificial_rewritten}
    \begin{aligned}
        \Pi_{-2} {f^{\varepsilon}}&= \sum_{k=0}^2 \Psi_k({\mathbf{v}})(\Psi_k, {f^{\varepsilon}}) \\
        \Pi_{-1} {f^{\varepsilon}}&= \sum_{k=0}^2 \Psi_k({\mathbf{v}})(\Psi_k, {f^{\varepsilon}}) + \Psi_1({\mathbf{v}})\Psi_2({\mathbf{v}})(\Psi_0, {f^{\varepsilon}}).
    \end{aligned}$$ We can thus view the linear kinetic models and used in [@LafitteMelisSamaey2017; @MelisSamaey2017], respectively, as a special simplified case of the linearized BGK equation.

Since the linearized BGK operator was shown to be nearly identical to the relaxation operators of the linear kinetic models in , we expect the spectral properties of the linearized BGK equation to closely resemble those in [@LafitteMelisSamaey2017] (for ${\nu}= 1$) or [@MelisSamaey2017] (for ${\nu}= \rho$). Therefore, it is expected that the construction of stable PI methods for the full BGK equation with ${\nu}= 1$ and stable TPI methods for with ${\nu}= \rho$ is practically identical to that in [@LafitteMelisSamaey2017] and [@MelisSamaey2017], respectively. The choice of method parameters for the full BGK equation will be discussed more closely in section \[sec:method\_param\].

Linearized Boltzmann equation and its spectrum {#subsec:lin_boltzmann_equation}
----------------------------------------------

### Linearization of the Boltzmann equation {#subsubsec:lin_boltzmann_equation}

To simplify the analysis of the Boltzmann equation , it is customary to linearize the collision operator ${{\mathcal{Q}}{\mathopen{}\mathclose\bgroup\originalleft}({f^{\varepsilon}}{\aftergroup\egroup\originalright})({\mathbf{v}})}$ around the global Maxwellian distribution ${{\mathcal{M}}^{{\rho^{\infty}},{{\mathbf{{\bar{v}}}}^{\infty}},{T^{\infty}}}_{{\mathbf{v}}}}= {\mathcal{M}}_{\mathbf{v}}^{1,0,1}$, which, for $D_v = 2$, is given by: $$\label{eq:global_maxwellian}
    {\mathcal{M}}_{\mathbf{v}}^{1,0,1} = \frac{1}{2\pi} \exp{\mathopen{}\mathclose\bgroup\originalleft}(-\frac{{{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}{\aftergroup\egroup\originalright}\vert}^2}{2}{\aftergroup\egroup\originalright}),$$ see, for instance, [@Cercignani1988; @Golse2005; @Saint-Raymond2009]. Subsequently, we consider small fluctuations of ${f^{\varepsilon}}$ around the global equilibrium , that is: $$\label{eq:linearization_f}
    {f^{\varepsilon}}({\mathbf{v}}) = {\mathcal{M}}_{\mathbf{v}}^{1,0,1}\big(1 + {g^{{\varepsilon}}}({\mathbf{v}})\big).$$ Using a similar shorthand notation as before, that is, $g' = {g^{{\varepsilon}}}({\mathbf{v}}')$, ${\mathcal{M}}' = {\mathcal{M}}_{{\mathbf{v}}'}^{1,0,1}$, and so on, we linearize the quadratic terms in as: $$\label{eq:linearizing_ff_ff}
    f'f_*' - ff_* = {\mathcal{M}}{\mathcal{M}}_*\big(g' + g_*' - g - g_*\big),$$ where we neglected second-order fluctuations in the second equality, and we used that ${\mathcal{M}}'{\mathcal{M}}_*' = {\mathcal{M}}{\mathcal{M}}_*$ for any Maxwellian distribution ${\mathcal{M}}$. By substituting and into the Boltzmann equation , and exploiting that ${\mathcal{M}}= {\mathcal{M}}_{\mathbf{v}}^{1,0,1}$ depends only on velocity ${\mathbf{v}}$, we obtain the *linearized Boltzmann equation* as: $$\label{eq:linearized_Boltzmann_equation}
    \partial_t {g^{{\varepsilon}}}+ {\mathbf{v}}\cdot \nabla_{{\mathbf{x}}} {g^{{\varepsilon}}}= \frac{1}{{\varepsilon}}{\mathcal{L}}{{g^{{\varepsilon}}}},$$ where ${\mathcal{L}}{{g^{{\varepsilon}}}}$ is the *linearized Boltzmann collision operator*, which reads: $$\label{eq:linearized_collision_operator}
    {\mathcal{L}}{{g^{{\varepsilon}}}} = \int_{{\mathbb{R}}^2}\int_{0}^{2\pi} B({{\mathopen{}\mathclose\bgroup\originalleft}\vert{{\mathbf{v}}_r}{\aftergroup\egroup\originalright}\vert},{\theta_\sigma}){\mathcal{M}}_*(g_*' + g' - g_* - g) d{\theta_\sigma}d{{\mathbf{v}}_{*}},$$ see also [@CIP; @Saint-Raymond2009]. Moreover, the operator ${\mathcal{L}}{{g^{{\varepsilon}}}}$ can be cast in the following form [@Cercignani1988 IV.5]: $$\label{eq:linearized_collision_operator_K_nu}
    {\mathcal{L}}{{g^{{\varepsilon}}}}({\mathbf{v}}) = {\mathcal{K}}{g^{{\varepsilon}}}({\mathbf{v}}) - {\nu}({{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}{\aftergroup\egroup\originalright}\vert}){g^{{\varepsilon}}}({\mathbf{v}}),$$ where ${\nu}({{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}{\aftergroup\egroup\originalright}\vert})$ is a local multiplication operator termed the *linearized collision frequency* that depends only on the magnitude of ${\mathbf{v}}$ and is defined as: $$\label{eq:collision_frequency_linearized}
    {\nu}({{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}{\aftergroup\egroup\originalright}\vert}) = \int_{{\mathbb{R}}^2}\int_{0}^{2\pi} B({{\mathopen{}\mathclose\bgroup\originalleft}\vert{{\mathbf{v}}_r}{\aftergroup\egroup\originalright}\vert},{\theta_\sigma}){\mathcal{M}}_* d{\theta_\sigma}d{{\mathbf{v}}_{*}},$$ and ${\mathcal{K}}$ is a non-local integral operator containing the remaining three terms in .

Using and writing ${f^{\varepsilon}}$ instead of ${g^{{\varepsilon}}}$, we rewrite as: $$\label{eq:linearized_Boltzmann_equation_framework}
    \partial_t {f^{\varepsilon}}+ {\mathbf{v}}\cdot \nabla_{\mathbf{x}}{f^{\varepsilon}}= -\frac{1}{{\varepsilon}}\big({\nu}({{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}{\aftergroup\egroup\originalright}\vert}){\mathcal{I}}- {\mathcal{K}}\big){f^{\varepsilon}}.$$ For inverse power law potentials under Grad’s cut-off assumption and for hard sphere particles, it can be proven that ${\mathcal{K}}$ is a compact operator on $L^2({\mathbb{R}}^2)$, see [@Cercignani1988; @Golse2005]. This implies that it maps the unit ball of ${\mathbb{R}}^2$ onto a finite-dimensional space [@Cercignani1988]. In that sense, it shares properties with the finite rank operators described previously. However, the linearized operator ${\mathcal{L}}= {\mathcal{K}}- {\nu}({{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}{\aftergroup\egroup\originalright}\vert}){\mathcal{I}}$ does not have finite rank. Furthermore, one could write formally: $${\mathcal{K}}= {\nu}({{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}{\aftergroup\egroup\originalright}\vert})(\Pi_\text{BGK} + {\mathcal{R}}),$$ for a certain remainder operator ${\mathcal{R}}$. Then, one can write as: $$\partial_t {f^{\varepsilon}}+ {\mathbf{v}}\cdot \nabla_{\mathbf{x}}{f^{\varepsilon}}= -\frac{{\nu}}{{\varepsilon}} ({\mathcal{I}}- \Pi_\text{BGK}){f^{\varepsilon}}- \frac{{\nu}}{{\varepsilon}} {\mathcal{R}}{f^{\varepsilon}},$$ which provides a connection with the linearized BGK projection operator.

### Linearized Boltzmann spectrum.

Analyzing the spectrum of the linearized Boltzmann collision operator ${\mathcal{L}}= {\mathcal{K}}- {\nu}({{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}{\aftergroup\egroup\originalright}\vert}){\mathcal{I}}$ in is more involved than in the BGK case. In general, the linearized Boltzmann collision operator has a spectrum that consists of (i) a non-empty essential (purely continuous) part that is entirely determined by the continuous spectrum of $-{\nu}({{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}{\aftergroup\egroup\originalright}\vert}){\mathcal{I}}$, and (ii) a set of discrete eigenvalues that is influenced by the operator ${\mathcal{K}}$, see, for instance, [@Cercignani1988] or [@BarangerMouhot2005; @Ellis1975]. In contrast, the spectrum of the linear kinetic relaxation operators and linearized BGK operator only consists of discrete eigenvalues. However, for Maxwellian particles with angular cut-off and, in particular, for pseudo-Maxwellian particles, it is known that the spectrum of ${\mathcal{L}}$ contains only discrete eigenvalues spread inside the interval $[-{\nu}(0),0]$ with ${\nu}({{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}{\aftergroup\egroup\originalright}\vert})$ given in [@Cercignani1988]. We have seen in Section \[subsubsec:lin\_boltzmann\_equation\] that the linearized version of the full Boltzmann equation reads: $$\partial_t {f^{\varepsilon}}+ {\mathbf{v}}\cdot \nabla_{{\mathbf{x}}} {f^{\varepsilon}}= -\frac{1}{{\varepsilon}}(\nu(|{\mathbf{v}}|){\mathcal{I}}- {\mathcal{K}}){f^{\varepsilon}}.$$ Let us apply the Fourier transform in the physical space: since the collision operator depends only on the velocity magnitude ${{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{v}}{\aftergroup\egroup\originalright}\vert}$, the only difference in the equation will be that the free transport term ${\mathbf{v}}\cdot \nabla_{\mathbf{x}}$ will become a multiplication operator (by $i \bm \gamma \cdot {\mathbf{v}}$, where $\bm \gamma$ is the spatial Fourier variable). One can then write the Fourier-transformed linear Boltzmann equation as: $$\label{eq:linearBoltzFourier}
    \partial_t {h^{{\varepsilon}}}= \frac 1{\varepsilon}{\mathcal{K}}{h^{{\varepsilon}}}- {\mathopen{}\mathclose\bgroup\originalleft}( \nu(|{\mathbf{v}}|)/{\varepsilon}+ i \, {\varepsilon}\bm \gamma \cdot {\mathbf{v}}{\aftergroup\egroup\originalright}) {h^{{\varepsilon}}},$$ where ${h^{{\varepsilon}}}$ is the Fourier transform in space of ${f^{\varepsilon}}$. Hence, the evolution of ${h^{{\varepsilon}}}$ is given by a compact perturbation of a (complex-valued) multiplication operator. It was proven in a series of papers that the spectrum of this Fourier-transformed collision operator has the following behavior as a function of $|\bm \gamma|$ and ${\varepsilon}$:

\[thm:linBoltzSpectrum\] The spectrum of the right hand side of equation consists of an essential part $\Sigma_e$ located to the left of a vertical line of negative real part and a discrete spectrum $\Sigma_d$ composed of:

-   *fast modes:* eigenvalues located at a distance at least $1/{\varepsilon}$ to the left of the imaginary axis;

-   *slow modes:* if $|{\varepsilon}| \ll 1$, there are *exactly* $4$ eigenvalues branches given by: $$\lambda^{(j)}(|\bm \gamma|) := i \, \lambda^{(j)}_1 \, {\varepsilon}|\bm \gamma| - \lambda^{(j)}_2 \, {\varepsilon}^2 |\bm \gamma|^2 + \mathcal O {\mathopen{}\mathclose\bgroup\originalleft}( {\varepsilon}^3 |\bm \gamma|^3 {\aftergroup\egroup\originalright}), \quad j \in \{0,1,2,3\},$$ for explicit constants $\lambda_1^{(j)} \in {\mathbb{R}}$ and $\lambda_2^{(j)}>0$.

A sketch of this result can be found in figure \[fig:spectrumBoltzFourier\].

![\[fig:spectrumBoltzFourier\] Spectrum of the Fourier transformed linearized Boltzmann operator, for small radial frequencies.](image/spectrum.pdf)

We observe in Theorem \[thm:linBoltzSpectrum\] that the discrete eigenvalues form a fast and a slow cluster, justifying the use of projective integration with easily computable parameters, in the spirit of [@LafitteMelisSamaey2017; @MelisSamaey2017]. Nevertheless, the presence of an essential spectrum is one of the reasons that one must use *telescopic* projective integration to solve the full Boltzmann equation: this spectral decomposition will give rise to new clusters of eigenvalues at the discrete level.

Note that one can mimic the proof of this result for the simpler BGK operator with $\nu = 1$ to obtain the same spectral behavior of the linearized operator without the essential part, justifying at the continuous level the results from [@LafitteMelisSamaey2017].

Method parameters for projective and telescopic projective integration {#sec:method_param}
----------------------------------------------------------------------

It still remains to select appropriate parameter values for the projective or telescopic projective integration methods. These are determined by ensuring that all eigenvalues of the kinetic problem under study fall within the stability region of the full projective method. In sections \[subsec:lin\_bgk\] and \[subsec:lin\_boltzmann\_equation\], we revealed that the spectra of the linearized kinetic equations either appear in two stationary eigenvalue disks (linearized BGK equation with ${\nu}= 1$) or are continuously spread along (a part of) the negative real axis (linearized BGK equation with ${\nu}= \rho$ and linearized Boltzmann equation). Since the construction of stable projective methods for both well-separated as well as continuously spread spectra is studied in previous works [@LafitteMelisSamaey2017] and [@MelisSamaey2017], respectively, we take over the main results here, which are summarized below.

### Stationary, well-separated spectrum

For the linearized BGK equation with $\nu = 1$, which falls into the class of kinetic models studied in [@LafitteMelisSamaey2017], it was shown in [@LafitteMelisSamaey2017] that the spectrum consists of two stationary, well-separated eigenvalue clusters (a fast and slow, dominant cluster). To accommodate these two clusters, the method parameters of projective integration can be selected such that its stability region splits up into two parts.

1.  First, the inner integrator time step ${\delta t}$ is chosen corresponding to the fastest time scale of the problem, which is of the order of ${\varepsilon}$. This centers one stability region of the projective method around the fast eigenvalues.

2.  Next, the number of inner integrator time steps $K$ is chosen such that all fast eigenvalues lie inside this stability region. In [@LafitteMelisSamaey2017], it was proven that we require $K \ge 2$.

3.  Last, the outer integrator time step ${\Delta t}$ is selected such that all dominant eigenvalues fall into the second stability region of the projective method.

Since both $K$ and ${\Delta t}$ are independent of the small-scale parameter ${\varepsilon}$, the resulting projective method has a cost that is also independent of ${\varepsilon}$, which becomes increasingly advantageous for ${{\varepsilon}\to 0}$.

### Continuously spread spectrum

When considering the linearized BGK equation with $\nu = \rho$, which was studied in [@MelisSamaey2017], the spectrum varies continuously over the negative real axis. This also holds true for the linearized Boltzmann equation, be it on a part of the negative real axis, see figure \[fig:spectrumBoltzFourier\]. In this case, we require that the stability region of the numerical method does not split up but instead comprises the entire negative real axis up to the fastest eigenvalue of the problem (a numerical method with this property is termed $\mathit{[0,1]}$*-stable*). Here, for simplicity, we assume that the fastest eigenvalue at $t = 0$ corresponds to the fastest possible eigenvalue for all other times $t > 0$. Since \[0,1\]-stable projective integration methods lose practically all of their potential speed-up, \[0,1\]-stable telescopic projective integration methods can be designed with much higher speed-ups.

1.  Similarly to projective integration, the innermost integrator time step $h_0$ of the telescopic projective integration method is chosen corresponding to the fastest time scale, which is of the order of ${\varepsilon}/\max_{x}\rho(x,0)$.

2.  Next, we fix the outermost time step we would like to use, taking into account a CFL-like stability constraint, as follows: $h_L = C{\Delta x}$.

3.  Before choosing the number of projective levels, we decide on the number of inner integrator time steps $K$, which we consider to be fixed on each projective level. For each chosen value of $K$ there is a corresponding maximal value of $M$ such that the stability region does not split up, see [@Gear2003telescopic] or [@MelisSamaey2017].

4.  The required number of projective levels $L$ to obtain a \[0,1\]-stable telescopic method is computed as (see [@MelisSamaey2017]): $$\label{eq:tpi_number_of_levels}
        L \approx \frac{\log(h_L) + \log(1/h_0)}{log(M+K+1)}.$$

5.  For the given values of $h_0$, $h_L$, $K$ and $L$ adapt the value of $M$ on the different projective levels such that the following equation: $$\label{eq:tpi_hL}
        h_L = \prod_{\ell=0}^{L-1} (M_\ell + K + 1)h_0$$ is valid (for further details, we refer the reader to [@MelisSamaey2017]).

For a \[0,1\]-stable telescopic method, the values of $h_L$, $M_\ell$ and $K$ are independent of ${\varepsilon}$. However, as indicated by equation , the number of projective levels increases as $O(\log(1/{\varepsilon}))$. As a consequence, the cost of a \[0,1\]-stable telescopic projective integration method is not completely ${\varepsilon}$-independent. However, the dependence is rather modest.

### Speedup {#subsubsec:speedup}

If we assume, as in [@Gear2003telescopic], that the overhead due to extrapolations is negligible and timestepping with the innermost integrator is computationally most demanding, the speedup ${\mathcal{S}}_L$ realized by the overall level-$L$ TPI method compared to naive forward Euler timestepping is given by: $$\label{eq:tpi_speedup} 
    {\mathcal{S}}_L = \prod_{\ell=0}^{L-1} \frac{M_\ell + K_\ell + 1}{K_\ell + 1},$$ that is, the ratio of the total number of naive forward Euler time steps within one outermost time step $h_L$ (see equation ) over the number of actual innermost steps in the TPI method.

Numerical experiments {#sec:results}
=====================

Here, we report simulation results for the BGK equation and the Boltzmann equation using pseudo-Maxwellian particles (that is, $\gamma = 0$ and $b_0$ constant). For each experiment, we shall compactly indicate the dimensions in space $(D_x)$ and velocity $(D_v)$ by writing “$D_x$D/$D_v$D" with $D_x,D_v \in \{1,2\}$. We begin with BGK in 1D/1D (Section \[subsec:bgk\_1d1d\]), and subsequently consider both the BGK and Boltzmann equation in 1D/2D (Section \[subsec:bgk\_boltzmann\_1d2d\]). Thereafter, we target a shock-bubble interaction problem for the BGK equation in 2D/2D (Section \[subsec:shock\_bubble\_interaction\]), and the more intricated Kelvin-Helmoltz-like instability in the same setting (Section \[subsec:KelvinHelmoltz\]). As a last experiment, we deal with the full Boltzmann equation in 2D/2D (Section \[subsec:boltzmann\_2d2d\]).

BGK in 1D/1D {#subsec:bgk_1d1d}
------------

As a first experiment, we focus on the nonlinear BGK equation in 1D/1D. We consider a Sod-like test case for $x \in [0,1]$ consisting of an initial centered Riemann problem with the following left and right state values: $$\label{eq:bgk_1d_sod}
    \begin{aligned}
        \begin{pmatrix} \rho_L \\ {\bar{v}}_L \\ T_L \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \\ 1 \end{pmatrix}, \qquad\qquad 
        \begin{pmatrix} \rho_R \\ {\bar{v}}_R \\ T_R \end{pmatrix} = \begin{pmatrix} 0.125 \\ 0 \\ 0.25 \end{pmatrix}.
    \end{aligned}$$ The initial distribution ${f^{\varepsilon}}(x,v,0)$ is then chosen as the Maxwellian corresponding to the above initial macroscopic variables. We impose outflow boundary conditions and perform simulations for $t \in [0,0.15]$. As velocity space, we take the interval $[-8,8]$, which we discretize on a uniform grid using $J=80$ velocity nodes. In all simulations, space is discretized using the WENO3 spatial discretization with ${\Delta x}= 0.01$. Below we regard three gas flow regimes, ${\varepsilon}= 10^{-1}$ (kinetic regime), ${\varepsilon}= 10^{-2}$ (transitional regime) and ${\varepsilon}= 10^{-5}$ (fluid regime), and for each regime, we compare solutions for two cases of collision frequency ${\nu}$ in the BGK equation , ${\nu}= 1$ and ${\nu}= \rho$.

#### Direct integration $\boldsymbol{({\varepsilon}= 10^{-1}}$ and $\boldsymbol{{\varepsilon}= 10^{-2})}$.

In the kinetic $({\varepsilon}= 10^{-1})$ and transitional $({\varepsilon}= 10^{-2})$ regimes, we compute the numerical solution for ${\nu}= 1$ and ${\nu}= \rho$ using the fourth order Runge-Kutta (RK4) time discretization with time step ${\delta t}= 0.1{\Delta x}$. The results are shown in figure \[fig:bgk\_1d1d\_nu\] for ${\nu}= 1$ (left) and ${\nu}= \rho$ (right), where we display the density $\rho$, macroscopic velocity ${\bar{v}}$ and temperature $T$ as given in at $t = 0.15$. In addition, we plot the heat flux $q$, which, in a general $D_v$-dimensional setting, is a vector ${\mathbf{q}}= {{\mathopen{}\mathclose\bgroup\originalleft}(q^d{\aftergroup\egroup\originalright})_{d=1}^{D_v}}$ with components given by: $$q^d = \frac{1}{2}\int_{{\mathbb{R}}^{D_v}} {{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{c}}{\aftergroup\egroup\originalright}\vert}^2c^d{f^{\varepsilon}}d{\mathbf{v}},$$ in which ${\mathbf{c}}= {{\mathopen{}\mathclose\bgroup\originalleft}(c^d{\aftergroup\egroup\originalright})_{d=1}^{D_v}} = {\mathbf{v}}- {\mathbf{{\bar{v}}}}$ is the peculiar velocity. The different regimes are shown by blue (kinetic) and purple (transitional) dots. The red line in each plot denotes the limiting $({{\varepsilon}\to 0})$ solution of each macroscopic variable, which all converge to the solution of the Euler system with ideal gas law $P = \rho T$ and heat flux $q = 0$.

[![\[fig:bgk\_1d1d\_nu\] Numerical solution of the BGK equation in 1D/1D with ${\nu}= 1$ (left) and ${\nu}= \rho$ (right) at $t = 0.15$ for a Sod-like shock test using the WENO3 scheme with ${\Delta x}= 0.01$. RK4 is used for ${\varepsilon}= 10^{-1}$ (blue dots) and ${\varepsilon}= 10^{-2}$ (purple dots) with ${\delta t}= 0.1{\Delta x}$. The PRK4 (left) and level-2 TPRK4 (right) methods are used for ${\varepsilon}= 10^{-5}$ (green dots). Red line: hydrodynamic limit solution $({{\varepsilon}\to 0})$. ](image/bgk1D1D_epsi125_nu1_nurho "fig:")]{}

#### Projective integration $\boldsymbol{({\varepsilon}= 10^{-5}}$ and $\boldsymbol{{\nu}= 1)}$.

In the fluid regime $({\varepsilon}= 10^{-5})$, direct integration schemes such as RK4 become too expensive due to a severe time step restriction, which is required to ensure stability of the method. Exploiting that the spectrum of the linearized BGK equation with ${\nu}= 1$ resembles that of the linear kinetic models used in [@LafitteMelisSamaey2017], see section \[sec:linearOpSpecProp\], we construct a projective integration method to accelerate time integration in the fluid regime. As inner integrator, we select the forward Euler time discretization with ${\delta t}= {\varepsilon}$. As outer integrator, we choose the fourth-order projective Runge-Kutta (PRK4) method, using $K=2$ inner steps and an outer step of size ${\Delta t}= 0.4{\Delta x}$. Figure \[fig:bgk\_1d1d\_nu\] (left) shows the macroscopic observables in the fluid regime for ${\nu}= 1$ at $t = 0.15$ (green dots). From this, we observe that the BGK solution is increasingly dissipative for increasing values of ${\varepsilon}$ since the rate with which ${f^{\varepsilon}}$ converges to its equilibrium ${{{\mathcal{M}}_{v}{\mathopen{}\mathclose\bgroup\originalleft}({f^{\varepsilon}}{\aftergroup\egroup\originalright})}}$ becomes slower. In contrast, for sufficiently small ${\varepsilon}$, relaxation to thermodynamic equilibrium occurs practically instantaneous and the Euler equations yield a valid description. Since this is a hyperbolic system, it allows for the development of sharp discontinuous and shock waves which are clearly seen in the numerical solution.

In this numerical test, the speed-up factor between a naive RK4 implementation and the projective integration method is $130.3$, namely formula with $L=1$ (1 projective level), $K_0=2$ and $M_0=397$.

[![\[fig:bgk\_rk4TPI\_comp\] Numerical solution of the BGK equation in 1D/1D with ${\nu}= \rho$, ${\varepsilon}= 10^{-5}$ at $t = 0.15$ for a Sod-like shock test using the WENO3 scheme with ${\Delta x}= 0.01$. Comparison between level-2 TPRK4 with ${\delta t}= 0.1{\Delta x}$ (solid blue line) and classical RK4 with a ${\Delta t}= 0.5{\varepsilon}$ (red dots).](image/bgk1D1D_epsi5_RK4TPIcomp "fig:")]{}

#### Telescopic projective integration $\boldsymbol{({\varepsilon}= 10^{-5}}$ and $\boldsymbol{{\nu}= \rho)}$.

Next, we repeat the above experiment taking ${\nu}= \rho$ in the BGK equation . We now design a $[0,1]$-stable telescopic projective integration method as in [@MelisSamaey2017] since, for this choice of ${\nu}$, the spectrum of the linearized BGK equation is spread along the negative real axis and is time-dependent. Therefore, the two-scale nature in case of ${\nu}= 1$ has become a multi-scale problem, which destroys the acceleration in time of projective integration. We construct a $[0,1]$-stable TPRK4 method consisting of $2$ projective levels with FE as innermost integrator with time step $h_0 = {\varepsilon}$, constant $K=6$ on each level and an outermost time step $h_2 = 0.4{\Delta x}$. The extrapolation step sizes $M$ on each level are calculated as $M=\{14.24, 11.83\}$. The results are shown by green dots in figure \[fig:bgk\_1d1d\_nu\] (right). We conclude that the effect of choosing ${\nu}= \rho$ primarily manifests itself in the transitional regime $({\varepsilon}= 10^{-2})$, for which the relaxation rate is not too slow nor too fast. Moreover, it is seen that this choice of collision frequency does not alter the hydrodynamic limit of the BGK equation, which is captured correctly by the telescopic scheme.

Finally, figure \[fig:bgk\_rk4TPI\_comp\] compares the solutions obtained with the level-2 TPRK4 method and a classical RK4 method with a very small time step (of order ${\varepsilon}$) for the stiff test case where ${\varepsilon}=10^{-5}$. We observe a very good agreement between the two simulations (only a very small difference can be seen in the heat flux), while the TPRK4 scheme is more than 10 times faster than the RK4 scheme, because of its bigger time steps. The former scheme would be even more efficient with smaller values of the relaxation parameter.

In this test, the speed-up factor between a naive RK4 implementation and the telescopic projective integration method is $8.2$, namely formula with $L=2$ (2 projective levels), $(K_0, K_1)=(6,6)$ and $(M_0,M1)=(14.24,11.83)$.

BGK and Boltzmann in 1D/2D {#subsec:bgk_boltzmann_1d2d}
--------------------------

The BGK equation was introduced as a simplified model for the Boltzmann equation capturing most essential features of the latter. Here, we investigate the difference between both models. Since the Boltzmann collision operator vanishes for a one-dimensional velocity space, in this section, we consider both models in 1D/2D. In the experiments, this is achieved by discretizing space $(x,y)$ on a grid of size $I_x \times 2$ and using homogeneous data along the $y$-direction such that the spatial derivative $\partial_y{f^{\varepsilon}}$ exactly cancels out.

We perform the Sod test of the previous section in 1D/2D, see also [@Filbet2010]. As velocity space, we take the domain $[-8,8]^2$, which we discretize on a uniform grid using $J_x = J_y = 32$ velocity nodes along each dimension. In all simulations, space is discretized using the WENO2 spatial discretization with ${\Delta x}= 0.01$. Below we regard two regimes, ${\varepsilon}= 10^{-2}$ (transitional regime) and ${\varepsilon}= 10^{-5}$ (fluid regime), and for each regime, we compare solutions for BGK with ${\nu}= 1$, BGK with ${\nu}= \rho$ and Boltzmann with pseudo-Maxwellian particles. To approximate the Boltzmann collision operator, we apply the fast spectral method described in section \[sub:FastSpectral\] using $N_\theta = 4$ discrete angles. This is enough because of the spectral accuracy of the trapezoidal rule applied to periodic functions (see [@FiMoPa2006] for more details on this topic).

#### Direct integration ($\boldsymbol{{\varepsilon}= 10^{-2}}$).

In the transitional regime, we perform all simulations using the RK4 method with time step ${\delta t}= 0.1{\Delta x}$, for which we display the results in figure \[fig:boltz\_bgk\_1d2d\_nu\] (left) for BGK with ${\nu}= 1$ (blue dots), BGK with ${\nu}= \rho$ (green dots) and the Boltzmann equation (red dots). From this, we observe that the BGK solution with ${\nu}= \rho$ is closer to the Boltzmann solution than the BGK solution with ${\nu}= 1$. This is as expected, since the BGK equation with ${\nu}= \rho$ correctly captures the loss term of the Boltzmann collision operator, see . Moreover, the discrepancy between Boltzmann and BGK with ${\nu}= \rho$ increases for higher order moments of ${f^{\varepsilon}}$; while the density (zeroth order moment) appears to coincide (to the naked eye), the heat flux (third order moment) reveals a clear difference between both models.

[![\[fig:boltz\_bgk\_1d2d\_nu\] Comparison between BGK and Boltzmann in 1D/2D for ${\varepsilon}= 10^{-2}$ (left) and ${\varepsilon}= 10^{-5}$ (right) at $t = 0.15$ for a Sod-like shock test using the WENO2 scheme with ${\Delta x}= 0.01$. Blue dots: BGK with ${\nu}= 1$; green dots: BGK with ${\nu}= \rho$; red dots: Boltzmann. For ${\varepsilon}= 10^{-2}$, we used RK4 with ${\delta t}= 0.1{\Delta x}$. For ${\varepsilon}= 10^{-5}$, we applied a PRK4 (blue dots) and a level-2 TPRK4 (green and red dots) method. ](image/boltz_bgk1D2D_epsi25_nu1_nurho "fig:")]{}

#### Projective methods ($\boldsymbol{{\varepsilon}= 10^{-5}}$).

In the fluid regime, the RK4 method becomes too expensive. To that end, for BGK with ${\nu}= 1$, we design a PRK4 method with FE as inner integrator using ${\delta t}= {\varepsilon}$, $K = 2$ inner steps and ${\Delta t}= 0.4{\Delta x}$. Due to the multi-scale nature of both the BGK relaxation operator with ${\nu}= \rho$ and the Boltzmann collision operator, we construct a $[0,1]$-stable level-2 TPRK4 method for both models using the FE scheme as innermost integrator with $h_0 = {\varepsilon}$. We set $K = 4$ constant on each level, compute the extrapolation step sizes as $M = \{14.24, 11.83\}$, and choose the outermost time step as $h_2 = 0.4{\Delta x}$. The results can be seen in figure \[fig:boltz\_bgk\_1d2d\_nu\] (right) using the same plotting style as in the left column. For all models, the projective and telescopic projective integration methods display the expected hydrodynamic limit.

In these numerical tests, the speed-up factor between a naive RK4 implementation and the telescopic projective integration method for the PRK4 method for the BGK model with constant relaxation is $133.3$, namely formula with $L=1$ (1 projective level), $K_0=2$ and $M_0=397$. The speedup for the TPRK4 method for the BGK model with nonconstant relaxation rate and the Boltzmann equation is $13$, namely formula with $L=2$ (2 projective levels), $(K_0, K_1)=(4,4)$ and $(M_0,M1)=(14.24,11.83)$.

Shock-bubble interaction in 2D/2D {#subsec:shock_bubble_interaction}
---------------------------------

Here, we consider the BGK equation in 2D/2D with constant collision frequency ${\nu}= 1$ and we investigate the interaction between a moving shock wave and a stationary smooth bubble, which was proposed in [@Torrilhon2006], see also [@Cai2010]. This problem consists of a shock wave positioned at $x = -1$ in a spatial domain ${\mathbf{x}}= (x,y) \in [-2,3] \times [-1,1]$ traveling with Mach number ${\mathit{Ma}}= 2$ into an equilibrium flow region. Over the shock wave, the following left $(x \le -1)$ and right $(x > -1)$ state values are imposed [@Cai2010]: $$\label{eq:bgk_2d_shock}
    \begin{aligned}
        \big(\rho_L, {\bar{v}}^x_L, {\bar{v}}^y_L, T_L\big) &= {\mathopen{}\mathclose\bgroup\originalleft}(\frac{16}{7}, \sqrt{\frac{5}{3}}\frac{7}{16}, 0, \frac{133}{64}{\aftergroup\egroup\originalright}) \\
        \big(\rho_R, {\bar{v}}^x_R, {\bar{v}}^y_R, T_R\big) &= (1, 0, 0, 1).
    \end{aligned}$$ Due to this initial profile, the shock wave will propagate rightwards into the flow region at rest $(x > -1)$. Moreover, in this equilibrium region, a smooth Gaussian density bubble centered at ${\mathbf{x}}_0 = (0.5,0)$ is placed, given by: $$\label{eq:bgk_2d_bubble}
    \rho({\mathbf{x}},0) = 1 + 1.5\exp{\mathopen{}\mathclose\bgroup\originalleft}(-16{{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{x}}- {\mathbf{x}}_0{\aftergroup\egroup\originalright}\vert}^2{\aftergroup\egroup\originalright}).$$ The initial density $\rho({\mathbf{x}},0)$ and temperature $T({\mathbf{x}},0)$ are visualized in figure \[fig:bgk\_2d2d\_shock\_bubble\_t00\]. Then, the initial distribution ${f^{\varepsilon}}({\mathbf{x}},{\mathbf{v}},0)$ is chosen as the Maxwellian corresponding to the initial macroscopic variables in -. We impose outflow and periodic boundary conditions along the $x$- and $y$-directions, respectively, and we perform simulations for $t \in [0,0.8]$. As velocity space, we take the domain $[-10,10]^2$, which we discretize on a uniform grid using $J_x = J_y = 30$ velocity nodes along each dimension. We discretize space using the WENO2 spatial discretization with $I_x = 200$ and $I_y = 25$. Furthermore, we consider a fluid regime by taking ${\varepsilon}= 10^{-5}$.

[![\[fig:bgk\_2d2d\_shock\_bubble\_t00\] Initial solution for density (top) and temperature (bottom) of the shock-bubble interaction in 2D/2D. The spatial domain $[-2,3] \times [-1,1]$ is discretized using $I_x = 200$ and $I_y = 25$. ](image/bgk2D2D_shock_bubble_t00 "fig:")]{}

Since we regard the BGK equation with constant collision frequency ${\nu}= 1$, the spectrum of the linearized BGK operator consists of two eigenvalue clusters. Therefore, we construct a projective integration method to speed up simulation in time. We select the PRK4 method with FE as inner integrator. The inner time step is fixed as ${\delta t}= {\varepsilon}$ and we use $K = 2$ inner steps in each outer integrator iteration. The outer time step is chosen as ${\Delta t}= 0.4{\Delta x}$. The simulated density and temperature at time $t = 0.8$ are displayed in figure \[fig:bgk\_2d2d\_shock\_bubble\_t08\]. We observe that the shock propagates in the positive $x$-direction and bumps into the stationary bubbly density. A similar behavior is seen for the temperature evolution.

To compare our results with those in [@Torrilhon2006], where the smallest value of ${\varepsilon}$ is chosen as ${\varepsilon}= 10^{-2}$, we regard the one-dimensional evolution of density and temperature along the axis $y = 0$. For $t \in \{0, 0.2, 0.4, 0.6, 0.8\}$, we plot these intersections in figure \[fig:bgk\_2d2d\_shock\_bubble\_y0\]. We conclude that we obtain the same solution structure at $t = 0.8$ as in [@Torrilhon2006]. However, our results are sharper and less dissipative supposedly due to the particular small value of ${\varepsilon}$ ($10^{-5}$ versus $10^{-2}$). In contrast to [@Cai2010], we nicely capture the swift changes in the temperature profile for $x \in [0.5,1]$ at $t = 0.8$.

Again, in this numerical test, the speed-up factor between a naive RK4 implementation and the projective integration method is $133.3$, namely formula with $L=1$ (1 projective level), $K_0=2$ and $M_0=397$.

[![\[fig:bgk\_2d2d\_shock\_bubble\_t08\] Numerical solution of the shock-bubble interaction at $t = 0.8$ using the BGK equation with ${\nu}= 1$ in 2D/2D. We discretized velocity space using $J_x = J_y = 30$. We applied a PRK4 method with FE as inner integrator and ${\delta t}= {\varepsilon}= 10^{-5}$ together with the WENO2 spatial discretization scheme with $I_x = 200$ and $I_y = 25$. ](image/bgk2D2D_shock_bubble_t08 "fig:")]{}

[![\[fig:bgk\_2d2d\_shock\_bubble\_y0\] Numerical solution of the shock-bubble interaction along $y = 0$ at $t = 0$ (black dashed), $t = 0.2$ (blue), $t = 0.4$ (purple), $t = 0.6$ (green) and $t = 0.8$ (red). ](image/bgk2D2D_shock_bubble_y0 "fig:")]{}

A Kelvin-Helmoltz like instability problem {#subsec:KelvinHelmoltz}
------------------------------------------

Keeping the same setting of the BGK equation in 2D/2D with constant collision frequency ${\nu}= 1$, we now turn to a less common test case in the field of collisional kinetic equation, the so-called Kelvin-Helmoltz instability. This phenomenon occurs when two fluids of different densities and in thermodynamic equilibrium move at different speeds. It is very well known that such a system will exhibit turbulent, unstable vortices at the interface between the two fluids, because of the velocity shear [@Helmoltz1868]. In order for these instabilities to develop, the Reynolds number of the fluids considered must be large. Using the von Karman relation [@sone2007molecular], which states that the Reynolds number is inversely proportional to the Knudsen number ${\varepsilon}$, we shall then choose a very small Knudsen number ${\varepsilon}= 5\cdot10^{-5}$ along with the following initial condition inspired from [@McNallyLyraPassy2012]: $$\label{eq:BGK2dKelvinHelmoltz}
    \begin{pmatrix} \rho_1 \\ {\bar{v}}^x_1 \\ {\bar{v}}^y_1 \\ T_1 \end{pmatrix} = \begin{pmatrix} 1 \\ 0.5 \\ 0.01 \sin{\mathopen{}\mathclose\bgroup\originalleft}(4 \pi x{\aftergroup\egroup\originalright}) \\ 1 \end{pmatrix} \quad (y \geq 0), \qquad\quad
    \begin{pmatrix} \rho_2 \\ {\bar{v}}^x_2 \\ {\bar{v}}^y_2 \\ T_2 \end{pmatrix} = \begin{pmatrix} 2 \\ -0.5 \\ 0.01 \sin{\mathopen{}\mathclose\bgroup\originalleft}(4 \pi x{\aftergroup\egroup\originalright})  \\ 1 \end{pmatrix} \quad (y < 0).$$ The initial distribution ${f^{\varepsilon}}({\mathbf{x}},{\mathbf{v}},0)$ is chosen as the Maxwellian corresponding to the initial macroscopic variables in . We impose periodic and outflow boundary conditions along the $x$- and $y$-directions, respectively, and we perform simulations for $t \in [0,1.6]$. As velocity space, we take the domain $[-8,8]^2$, which we discretize on a uniform grid using $J_x = J_y = 30$ velocity nodes along each dimension. We discretize space using the WENO2 method on $[-0.5,0.5]\times[-0.5,0.5]$ with $I_x = I_y= 100$.

Since we consider again the BGK equation with constant collision frequency ${\nu}= 1$, the spectrum of the linearized BGK operator consists of two eigenvalue clusters. Therefore, we construct a projective integration method to speed up simulation in time. We select the PRK4 method with FE as inner integrator. The inner time step is fixed as ${\delta t}= {\varepsilon}$ and we use $K = 3$ inner steps in each outer integrator iteration (this test case is stiffer than the previous one because of the turbulent regime). The outer time step is chosen as ${\Delta t}= 0.45{\Delta x}$. The simulated density and pressure at time $t = 0.4$, $0.9$ and $1.6$ are displayed in figure \[fig:bgk\_2d2d\_KH\_Density-Pressure\], along with the vector field $ {\mathbf{{\bar{v}}}}$ with the contour lines of density at time $t=0.9$ in figure \[fig:bgk\_2d2d\_KH\_velocityField\].

In this numerical test, the speed-up factor between a naive RK4 implementation and the projective integration method is $22.5$, namely formula with $L=1$ (1 projective level), $K_0=3$ and $M_0=86$.

[![\[fig:bgk\_2d2d\_KH\_Density-Pressure\] Density (left) and pressure (right) of the Kelvin-Helmoltzm-like instability at times $t=0.6$ (first row), $0.9$ (second row) and $1.6$ (third row). ](image/bgk2D2D_KH_Density-Pressure "fig:")]{}

[![\[fig:bgk\_2d2d\_KH\_velocityField\]Velocity field and density line of the Kelvin-Helmoltzm-like instability at time $t=0.9$.](image/bgk2D2D_KH_velocityField "fig:")]{}

As expected, we observe vortices developing along the velocity shear line $y=0$, and then expanding, forming small scale structures as time evolves, hence validating the turbulent behavior exhibited by the fluid for this Reynolds number. Note that the use of high-order methods in both space (WENO2) and time (PRK4) is crucial, because of these fine scale structures. Choosing less accurate (and more dissipative) schemes, such as forward Euler in time and upwind in space, leads to a uniform density and pressure instead of vortices in our numerical tests.

Boltzmann in 2D/2D {#subsec:boltzmann_2d2d}
------------------

As a last experiment, we concentrate on the Boltzmann equation with pseudo-Maxwellian particles in 2D/2D. As initial configuration for the gas, we consider the double Sod shock test, that is, for ${{\mathbf{x}}= (x,y) \in [-0.5,0.5]^2}$, we set: $$\label{eq:boltzmann_2d_sod}
    \begin{aligned}
        \begin{pmatrix} \rho_1 \\ {\mathbf{{\bar{v}}}}_1 \\ T_1 \end{pmatrix} = \begin{pmatrix} 0.1 \\ \boldsymbol{0} \\ 1 \end{pmatrix} \quad (xy \le 0), \qquad\quad
        \begin{pmatrix} \rho_2 \\ {\mathbf{{\bar{v}}}}_2 \\ T_2 \end{pmatrix} = \begin{pmatrix} 1 \\ \boldsymbol{0} \\ 1 \end{pmatrix} \quad (\text{otherwise}).
    \end{aligned}$$ The initial distribution ${f^{\varepsilon}}({\mathbf{x}},{\mathbf{v}},0)$ is then chosen as the Maxwellian corresponding to the above macroscopic variables. We impose outflow boundary conditions along both dimensions and perform simulations for $t \in [0,0.16]$. As velocity space, we take the domain $[-8,8]^2$, which we discretize on a uniform grid with $J_x = J_y = 32$ velocity nodes along each dimension. Furthermore, we discretize the spatial domain using the WENO2 spatial discretization with $I_x = I_y = 64$ grid points along each dimension, and we fix ${\varepsilon}= 5 \cdot 10^{-5}$ (that is, we are in the fluid regime). The Boltzmann collision operator is approximated again using the fast spectral method described in section \[sub:FastSpectral\] with $N_\theta = 4$ discrete angles.

For the time simulation of the Boltzmann equation, we apply a level-2 TPRK4 method with FE as innermost integrator using $h_0 = {\varepsilon}$ as innermost time step. We set $K = 3$ constant on each level and compute the extrapolation step sizes as $M = \{6.66, 4.80\}$. The outermost time step is chosen as ${\Delta t}= 0.3{\Delta x}$. In figure \[fig:boltz\_2d2d\_sod\], we plot various macroscopic observables of interest at $t = 0.16$. The density, macroscopic velocity along $x$ and temperature are computed as the moments of ${f^{\varepsilon}}$, see . Then, pressure, energy and the Mach number are obtained, respectively, as: $$P = \rho T, \qquad\quad E = \frac{1}{2}\rho{{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{{\bar{v}}}}{\aftergroup\egroup\originalright}\vert}^2 + P, \qquad\quad {\mathit{Ma}}= \frac{{{\mathopen{}\mathclose\bgroup\originalleft}\vert{\mathbf{{\bar{v}}}}{\aftergroup\egroup\originalright}\vert}}{\sqrt{T}}.$$ In this test, the speed-up factor between a naive RK4 implementation and the telescopic projective integration method is $5.9$, namely formula with $L=2$ (2 projective levels), $(K_0, K_1)=(3,3)$ and $(M_0,M1)=(6.66,4.80=)$.

[![\[fig:boltz\_2d2d\_sod\] Numerical solution of the Boltzmann equation with pseudo-Maxwellian particles in 2D/2D at $t = 0.16$ for a double Sod shock test . Velocity space is discretized using $J_x = J_y = 32$. We applied a level-2 TPRK4 method with FE as innermost integrator and $h_0 = {\varepsilon}= 5 \cdot 10^{-5}$ together with WENO2 with $I_x = I_y = 64$. ](image/boltz_2D2D_sod "fig:")]{}

Conclusions {#sec:conclusions}
===========

We extended projective and telescopic projective integration methods, studied in [@LafitteMelisSamaey2017; @MelisSamaey2017] for kinetic BGK-like equations, to allow for a fully explicit, high-order time simulation of both the nonlinear BGK equation and the full Boltzmann equation for pseudo-Maxwellian particles. We developed a framework of linearized operators, which revealed that the linearized BGK operator closely resembles the relaxation operator of the linear kinetic models used in [@LafitteMelisSamaey2017; @MelisSamaey2017]. As a result, the design of stable projective and telescopic projective integration methods for simulation of the nonlinear BGK equation is practically identical to that presented in [@LafitteMelisSamaey2017] and [@MelisSamaey2017], respectively. Since the spectrum of the linearized Boltzmann equation for pseudo-Maxwellian particles is known to be spread along the negative real axis, accelerated explicit time integration of the full Boltzmann equation required the construbtion of $[0,1]$-stable TPI methods. Although such methods are not completely asymptotic-preserving, the cost scales only logarithmically with the stiffness of the problem, which we consider acceptable. We designed and applied the projective methods to a number of model problems of increasing dimension and complexity, showing the potential of these schemes, with speed-ups that ranges from 8 for the most complicated Boltzmann model to more than a hundred for the BGK equation with linear relaxation rate.

Because of their explicit nature and hierarchical structure, telescopic projective integration methods naturally lend themselves to adaptivity. In future work, we will explore the use of space-dependent hierarchies of projective levels, adapted to the local values of ${\varepsilon}$, to further increase computational efficiency.

Free transport with WENO scheme {#sub:weno}
===============================

In weighted essentially non-oscillatory (WENO) methods, introduced in [@Liu1994], the stencil is built adaptively for each finite volume cell ${\mathcal{C}}_i = [x_{i-1/2},x_{i+1/2}]$. A WENO method considers all $k$ possible stencils for cell ${\mathcal{C}}_i$ at once and computes a weighted average over the $k$ resulting reconstructed solutions. The $k$ stencils for cell ${\mathcal{C}}_i$ are given by: $$\label{eq:weno_stencils}
    {\mathcal{S}}_i^\ell = \{{\mathcal{C}}_{i-\ell}, \ldots, {\mathcal{C}}_{i-\ell+k-1}\} \qquad (\ell = 0, \ldots, k-1),$$ in which the left shift $\ell$ is used to iterate over the possible stencils. As in [@Shu1998], for each stencil in , the reconstructed solution on the left side of $x_{i+1/2}$ and on the right side of $x_{i-1/2}$ are calculated as: $$\label{eq:weno_reconstructed_solutions}
    \begin{aligned}
        u_{i+1/2}^\ell &= \sum_{r=0}^{k-1} c_{\ell,r}U_{i-\ell+r} \\
        u_{i-1/2}^\ell &= \sum_{r=0}^{k-1} c_{\ell-1,r}U_{i-\ell+r}
    \end{aligned}
    \qquad\quad (\ell = 0, \ldots, k-1),$$ where we dropped the minus and plus superscripts in , respectively, for clarity. In what follows, we concentrate on the right boundary of ${\mathcal{C}}_i$. The left boundary is treated analogously.

Using the $k$ reconstructed solutions in , we could define the following weighted average: $$\label{eq:weno_reconstruction_averaged_smooth}
    u_{i+1/2}^- = \sum_{\ell=0}^{k-1} d_{\ell}u_{i+1/2}^\ell,$$ for some weights $d_{\ell}$, $\ell=0, \ldots, k-1$. Equation can be seen as the reconstruction of the solution on the left side of $x_{i+1/2}$ using a stencil of $2k-1$ consecutive cells. The weights $d_{\ell}$ can be chosen such that has an order of accuracy of $2k-1$, that is: $$u_{i+1/2}^- = u(x_{i+1/2}) + O{\mathopen{}\mathclose\bgroup\originalleft}({\Delta x}^{2k-1}{\aftergroup\egroup\originalright}),$$ which holds if the solution is smooth over all cells of . In [@Arandiga2011] it is proven that the weights $d_{\ell}$ can be written in terms of binomial coefficients as: $$d_{\ell} = \dfrac{\dbinom{k-1}{k-1-\ell}\dbinom{k}{k-1-\ell}}{\dbinom{2k-1}{k}} \qquad (\ell = 0, \ldots, k-1),$$ from which it is straightforward to see that $0 < d_{\ell} \le 1$ and $\sum_{\ell=0}^{k-1} d_{\ell} = 1$ such that represents a convex combination of the $k$ reconstructed solutions. For instance, we compute: $$\label{eq:weno_dl}
    \begin{alignedat}{4}
        k &=1: &&\quad d_0 = 1&& \\
        k &=2: &&\quad d_0 = \frac{2}{3}, &&\quad d_1 = \frac{1}{3} \\
        k &=3: &&\quad d_0 = \frac{3}{10}, &&\quad d_1 = \frac{6}{10}, \quad d_2 = \frac{1}{10}.
    \end{alignedat}$$ The idea of WENO methods is to use the weights $d_{\ell}$ only for stencils over which the solution is smooth and allocate very small weights to stencils containing one or more discontinuities. If we denote these WENO weights by $\omega_{\ell}$, we write the WENO reconstruction as: $$    u_{i+1/2}^- = \sum_{\ell=0}^{k-1} \omega_{\ell}u_{i+1/2}^\ell,$$ for which we require that the weights $\omega_{\ell}$ also form a convex set, that is: $$\label{eq:weno_omegal_convex}
    0 < \omega_\ell \le 1, \qquad\quad
    \sum_{\ell=1}^{k-1}\omega_\ell = 1.$$ In [@Liu1994] the second constraint in is guaranteed by writing $\omega_\ell$ in terms of constants $\alpha_\ell$ as: $$    \omega_\ell = \frac{\alpha_\ell}{\displaystyle\sum_{s=0}^{k-1} \alpha_s} \qquad (\ell = 0, \ldots, k-1).$$ The constants $\alpha_\ell$ are expressed in terms of the weights $d_{\ell}$, as given in , and new constants $\beta_{\ell}$, leading to: $$\label{eq:weno_alphal}
    \alpha_\ell = \frac{d_\ell}{(\delta + \beta_\ell)^2} \qquad (\ell = 0, \ldots, k-1),$$ in which the coefficients $\beta_\ell$ determine the smoothness of each stencil and are therefore termed the *smoothness indicators*. Since the smoothness indicators can become zero, a small constant $\delta > 0$ is added in that avoids division by zero. A typical choice is $\delta = 10^{-6}$ and the numerical experiments in [@Jiang1996] suggest that the resulting WENO method is not sensitive to the choice of $\delta$ as soon as $10^{-7} \le \delta \le 10^{-5}$.

The critical step in the design of WENO methods appears to be the definition of the smoothness indicators $\beta_{\ell}$. In [@Jiang1996], the authors proposed the following smoothness indicators: $$\label{eq:weno_smoothness_indicators}
    \beta_\ell = \sum_{s=1}^{k-1}{\Delta x}^{2s-1} \int_{x_{i-1/2}}^{x_{i+1/2}} {\mathopen{}\mathclose\bgroup\originalleft}(\frac{d^s p_i^\ell(x)}{dx^s}{\aftergroup\egroup\originalright})^2 dx \qquad (\ell = 0, \ldots, k-1).$$ Each smoothness indicator in includes the information of all derivatives that are approximated by the reconstruction polynomial $p_i^\ell(x)$ based on stencil ${\mathcal{S}}_i^\ell$ as a measure of smoothness of the true solution in cell ${\mathcal{C}}_i$. The factor ${\Delta x}^{2s-1}$ in front of the integral ensures that $\beta_\ell$ is independent of ${\Delta x}$ when computing the derivative of $p_i^\ell(x)$ in . For instance, for $k = 2$, the smoothness indicators in are calculated as: $$\label{eq:weno_betal_k2}
    \begin{aligned}
        \beta_0 &= (U_{i} - U_{i+1})^2 \\
        \beta_1 &= (U_{i-1} - U_{i})^2.
    \end{aligned}$$ For higher-order reconstructions, the derivatives in can be computed by using a general expression of the reconstruction polynomial. For $k=3$, after rearranging terms, we find (see [@Jiang1996; @Shu1998]): $$\label{eq:weno_betal_k3}
    \begin{aligned}
        \beta_0 &= \frac{1}{4}(3U_{i} - 4U_{i+1} + U_{i+2})^2 + \frac{13}{12}(U_{i} - 2U_{i+1} + U_{i+2})^2 \\
        \beta_1 &= \frac{1}{4}(U_{i-1} - U_{i+1})^2 + \frac{13}{12}(U_{i-1} - 2U_{i} + U_{i+1})^2 \\
        \beta_2 &= \frac{1}{4}(U_{i-2} - 4U_{i-1} + 3U_{i})^2 + \frac{13}{12}(U_{i-2} - 2U_{i-1} + U_{i})^2,
    \end{aligned}$$ see also the work of [@Henrick2005]. For orders $k=4,5,6$, the corresponding expressions were obtained in [@Balsara2000] and extended to even higher orders $k=7,8,9$ in [@Gerolymos2009]. From the expressions in and , we see that $\beta_\ell = O({\Delta x}^2)$ as soon as the solution is smooth over the stencil, while we have $\beta_\ell = O(1)$ for stencils containing a discontinuity. Using equation , this results in weights $\omega_\ell \approx d_\ell$ for smooth solutions and $\omega_\ell = O({\Delta x}^4)$ for discontinuous solutions, as desired.

Evaluating the Boltzmann collision operator with a fast spectral scheme {#sub:FastSpectral}
=======================================================================

The fast spectral discretization of the Boltzmann operator taken from [@MoPa:2006] employed in this work is described in this appendix.

To this aim, and since the Boltzmann collision operator acts only on the velocity variables, we focus on a given spatial cell ${\mathbf{x}}_j$ at a given instant of time $t^n$. Hence, only the dependency on the velocity variable ${\mathbf{v}}$ is considered for the distribution function $f$, i.e. $f=f({\mathbf{v}})$.

The first step to construct our spectral discretization is to truncate the integration domain of the Boltzmann integral . As a consequence, we suppose the distribution function $f$ to have compact support on the ball $\Ball_0(R)$ of radius $R$ centered in the origin. Since one can prove (see e.g. [@PaRuSINUM2000]) that $\supp (Q(f)(v)) \subset \Ball_0({\sqrt 2}R)$, in order to write a spectral approximation which avoids aliasing, it is sufficient that the distribution function $f({\mathbf{v}})$ is restricted on the cube $[-T,T]^{D_v}$ with $T \geq (2+{\sqrt 2})R$. Successively, one should assume $f({\mathbf{v}})=0$ on $[-T,T]^{D_v} \setminus \Ball_0(R)$ and extend $f$ to a periodic function on the set $[-T,T]^{D_v}$. Let observe that the lower bound for $T$ can be improved. For instance, the choice $T=(3+{\sqrt 2})R/2$ guarantees the absence of intersection between periods where $f$ is different from zero. However, since in practice the support of $f$ increases with time, we can just minimize the errors due to aliasing [@canuto:88] with spectral accuracy.

To further simplify the notation, let us take $T=\pi$ and hence $R=\lambda\pi$ with $\lambda = 2/(3+\sqrt{2})$ in the following. We denote by $\QL_B(f)$ the Boltzmann operator with cut-off. Hereafter, using one (bold) index to denote the $D_v$-dimensional sums, we have that the approximate function $f_N$ can be represented as the truncated Fourier series by $$f_N({\mathbf{v}}) = \sum_{{\mathbf{k}}=-N/2}^{N/2} \f_{\mathbf{k}}e^{i {\mathbf{k}}\cdot {\mathbf{v}}},
\label{eq:FU}$$ where the ${\mathbf{k}}^{th}$ Fourier coefficient is given by $$\f_{\mathbf{k}}= \frac{1}{(2\pi)^{D_v}}\int_{[-\pi,\pi]^{D_v}} f({\mathbf{v}})
e^{-i {\mathbf{k}}\cdot {\mathbf{v}}}\,d{\mathbf{v}}.$$ We then obtain a spectral quadrature of our collision operator by projecting on the space of trigonometric polynomials of degree less or equal to $N$, i.e. $${\hat {{\mathcal{Q}}}}_{\mathbf{k}}=\int_{[-\pi,\pi]^{D_v}}
\QL_B(f_N) \,
e^{-i {\mathbf{k}}\cdot {\mathbf{v}}}\,dv, \quad {\mathbf{k}}=-N/2,\ldots,N/2. 
\label{eq:VAR}$$ Finally, by substituting expression in one gets after some computations $${\hat {{\mathcal{Q}}}}_{\mathbf{k}}= \sum_{\substack{{\mathbf{l}},{\mathbf{m}}=-N/2\\ {\mathbf{l}}+{\mathbf{m}}={\mathbf{k}}}}^{N/2} \f_{\mathbf{l}}\,\f_{\mathbf{m}}\,
\bb({\mathbf{l}},{\mathbf{m}}),\quad {\mathbf{k}}=-N,\ldots,N,
\label{eq:CF1}$$ where $\bb({\mathbf{l}},{\mathbf{m}})=\B({\mathbf{l}},{\mathbf{m}})-\B({\mathbf{m}},{\mathbf{m}})$ are given by $$\B({\mathbf{l}},{\mathbf{m}}) = \int_{\Ball_0(2\lambda\pi)}\int_{\mathbb{S}^{D_v-1}} 
B(|q|, \cos\theta) e^{-i({\mathbf{l}}\cdot q^++{\mathbf{m}}\cdot q^-)}\,d\omega\,dq.$$ with $$q^{+} = \frac12(q+\vert q\vert \omega), \quad
q^{-} = \frac12(q-\vert q\vert \omega).$$ Let us notice that the naive evaluation of (\[eq:CF1\]) requires $O(n^2)$ operations, where $n=N^3$. This causes the spectral method to be computationally very expensive, especially in dimension three. In order to reduce the number of operations needed to evaluate the collision integral, the main idea is to use another representation of , the so-called Carleman representation [@Carl:EB:32] which is obtained by using the following identity $$\frac{1}{2} \, \int_{\mathbb{S}^{D_v-1}} F(|u|\sigma - u) \, d\sigma     = \frac{1}{|u|^{d-2}} \, \int_{\mathbb{R}^{D_v}} \delta(2 \, x \cdot u + |x|^2) \, F(x) \, dx.$$ This gives in our context for the Boltzmann integral $${{\mathcal{Q}}}(f)= \int_{{\mathbb{R}}^{D_v}} \int_{{\mathbb{R}}^{D_v}} {\tilde B}(x,y) 
 \delta(x \cdot y) 
 {\mathopen{}\mathclose\bgroup\originalleft}[ f(v + y) \, f(v+ x) - f(v+x+y) \, f(v) {\aftergroup\egroup\originalright}] \, dx \,
 dy,$$ with $$\label{eq:Btilde}
  \tilde{B}(|x|,|y|) =
  2^{D_v-1} \, \sigma{\mathopen{}\mathclose\bgroup\originalleft}(\sqrt{|x|^2+|y|^2}, \frac{|x|}{\sqrt{|x|^2+|y|^2}} {\aftergroup\egroup\originalright}) \, (|x|^2+|y|^2)^{-\frac{D_v-2}2}.$$ This transformation yields the following new spectral quadrature formula $$\label{eq:ode}
 \hat{{{\mathcal{Q}}}}_{\mathbf{k}}=
 \sum_{\underset{{\mathbf{l}}+{\mathbf{m}}={\mathbf{k}}}{{\mathbf{l}},{\mathbf{m}}=-N/2}}^{N/2} {\hat{\beta}}_F({\mathbf{l}},{\mathbf{m}}) \, \hat{f}_{\mathbf{l}}\, \hat{f}_{\mathbf{m}}, \ \ \
 {\mathbf{k}}=-N,...,N$$ where ${\hat{\beta}}_F({\mathbf{l}},{\mathbf{m}})=\B_F({\mathbf{l}},{\mathbf{m}})-\B_F({\mathbf{m}},{\mathbf{m}})$ are now given by $$\B_F({\mathbf{l}},{\mathbf{m}}) = \int_{\Ball_0(R)} \int_{\Ball_0(R)}
 \tilde{B}(x,y) \, \delta(x \cdot y) \, 
 e^{i ({\mathbf{l}}\cdot x+ {\mathbf{m}}\cdot y)} \, dx \, dy.$$ Now, in order to reduce the number of operation needed to evaluate , we look for a convolution structure. The aim is to approximate each ${\hat{\beta}}_F({\mathbf{l}},{\mathbf{m}})$ by a sum $${\hat{\beta}}_F({\mathbf{l}},{\mathbf{m}}) \simeq \sum_{p=1} ^{A} \alpha_p ({\mathbf{l}}) \alpha' _p ({\mathbf{m}}),$$ where $A$ represents the number of finite possible directions of collisions. This finally gives a sum of $A$ discrete convolutions and, consequently, the algorithm can be computed in $O(A \, N \log_2 N)$ operations by means of standard FFT technique [@canuto:88].

In order to get this convolution form, we make the decoupling assumption $$\tilde{B}(x,y) = a(|x|) \, b(|y|).$$ This assumption is satisfied if $\tilde B$ is constant. This is the case of Maxwellian molecules in dimension two, which is the case we shall consider for the numerical simulations of section \[sec:results\]. Indeed, using the kernel in , one has $$\tilde{B}(x,y) = 2^{D_v - 1} b_0 (|x|^2+|y|^2)^{-\frac{D_v-\alpha-2}2},$$ so that $\tilde{B}$ is constant if $D_v = 2$ and $\alpha = 0$. Here we write $x$ and $y$ in spherical coordinates $x = \rho e$ and $y = \rho' e'$ to get $$\B_F({\mathbf{l}},{\mathbf{m}}) = \frac14 \, \int_{\mathbb{S}^1} \int_{\mathbb{S}^1}
 \delta(e \cdot e') \,
 {\mathopen{}\mathclose\bgroup\originalleft}[ \int_{-R} ^R e^{i \rho ({\mathbf{l}}\cdot e)} \, d\rho {\aftergroup\egroup\originalright}] \,
 {\mathopen{}\mathclose\bgroup\originalleft}[ \int_{-R} ^R e^{i \rho' ({\mathbf{m}}\cdot e')} \, d\rho' {\aftergroup\egroup\originalright}] \, de \,
 de'.$$ Then, denoting $ \phi_R ^2 (s) = \int_{-R} ^R e^{i \rho s} \, d\rho,$ for $s \in {\mathbb{R}}$, we have the explicit formula $$\phi_R ^2 (s) = 2 \, R \,{\operatorname{Sinc}} (R s),$$ where ${\operatorname{Sinc}}(x)=\frac{\sin(x)}{x}$. This explicit formula is further plugged in the expression of $\B_F(l,m)$ and using its parity property, this yields $$\B_F ({\mathbf{l}},{\mathbf{m}}) =  \int_0 ^{\pi} \phi_R ^2 ({\mathbf{l}}\cdot e_{\theta})\, \phi_R ^2 ({\mathbf{m}}\cdot e_{\theta+\pi/2}) \,
 d\theta.$$ Finally, a regular discretization of $N_\theta$ equally spaced points, which is spectrally accurate because of the periodicity of the function [@kurganov2007spectral], gives $$\B_F ({\mathbf{l}},{\mathbf{m}}) = \frac{\pi}{M} \, \sum_{p=1} ^{N_\theta} \alpha_p ({\mathbf{l}}) \alpha' _p ({\mathbf{m}}),$$ with $$\alpha _p ({\mathbf{l}}) = \phi_R ^2 ({\mathbf{l}}\cdot e_{\theta_p}), \hspace{0.8cm} \alpha' _p ({\mathbf{m}}) = \phi_R ^2 ({\mathbf{m}}\cdot e_{\theta_p+\pi/2})$$ where $\theta_p = \pi p/N_\theta$.

[^1]: Department of Computer Science, K.U. Leuven, Celestijnenlaan 200A, 3001 Leuven, Belgium ([ward.melis@cs.kuleuven.be]{}).

[^2]: Laboratoire Paul Painlevé, Université de Lille, Cité Scientifique, 59655 Villeneuve d’Ascq, France([thomas.rey@math.univ-lille1.fr]{})

[^3]: Department of Computer Science, K.U. Leuven, Celestijnenlaan 200A, 3001 Leuven, Belgium ([giovanni.samaey@cs.kuleuven.be]{}).
