Keywords: space-homogeneous Boltzmann equation; tensor train; cross approximation; conservation

1 Introduction↩︎

The Boltzmann equation is a fundamental equation in kinetic theory that describes the evolution of the probability density function (pdf) of a dilute gas through binary collisions. Its central computational challenge lies in the fast computation of non-linear and non-local collision operator, a five-dimensional integral that is notoriously expensive to evaluate. There have been two classes of numerical methods for solving the Boltzmann equation: deterministic and stochastic. Stochastic methods are easy to implement but suffer from slow convergence and noise. In this paper, we will focus on the deterministic methods. Among the existing works, the fast Fourier spectral method [4]–[7] is popular because it can reduce the computational cost to \(\mathcal{O}(Mn^4 \log n)\) [7] by exploiting the convolution structure of the collision operator, where \(M\) is the number of angular quadrature points and \(n\) is the number of grid points in each velocity direction. In [8] the authors proposed a method to enforce conservation of mass, momentum and energy based on Lagrange multiplier. The Petrov-Galerkin approach expresses the solution as combinations of spherical harmonics and Hermite/Laguerre polynomials [9]–[11]. The solution is defined on the whole \(\mathbb{R}^{3}\) space, hence there is no need to choose a cut-off domain. In [12], the authors proposed a new approach that significantly reduces the memory cost and time cost, by taking advantage of the symmetry of collision operator. A different line of work is represented by discrete velocity methods [13]–[15], which discretize the velocity variable directly and are valued for their robustness and their ability to preserve structural properties such as conservation, positivity, and entropy dissipation, at the price of higher computational cost and typically lower-order accuracy. However, the aforementioned work does not provide mechanism to deal with the curse of dimensionality inherent in the high-dimensional phase space of the Boltzmann equation.

To tackle the high dimensionality of the Boltzmann equation, in recent years, there have been tremendous development of low-rank methods, which explores the low-rank properties of the probability density function [16]. We mention [17]–[20] which computes low-rank solutions to Boltzmann or BGK equations, among many others. However, most existing work focuses on the low-rank property between the velocity and the spatial variables, and the low-rank structure in the velocity space has not been explored for the Boltzmann collision operator due to the complex structure of the collision operator. On the other hand, we note the low-rank property is well motivated because the equilibrium solution is low-rank in velocity space.

In this paper, we aim to develop a fast method to compute the collisional integral term exploring the low-rank property of the solution in the velocity space. The work is motivated by the recent seminal contribution in the analysis of Boltzmann and Landau equations [1]–[3], where a lifting-projection framework is introduced for analyzing the well-posedness of those equations.

We present a numerical version of this procedure: the lifting-projection (LP) scheme for the Boltzmann equation. The key idea is to lift the original 3D equation to the 6D linear Kac master equation, evolve for one time step, and project the solution back to marginal distribution. The second contribution is a low-rank tensor method for evaluating the collision operator which computes with the lifted solution in 6D. We use the tensor train (TT) format [21] and present a numerical scheme in TT format utilizing interpolation and cross approximation, enhanced with conservation correction. The total computational cost per time step is \(\mathcal{O}(nMmR^2 r^2 X + nR^3X)\), where \(m\) is the width of interpolation stencil which is associated to the order of accuracy of the scheme, \(M\) is the number of spherical quadrature points, \(R\) and \(r\) are the max TT ranks of the lifted and original pdf, and \(X\) is the number of sweeps in cross approximation. We observe numerically that \(X\approx R\). Therefore, we can see that the computational cost scales linearly with \(n,\) and when the ranks \(R, r\) are small, this represents a dramatic computational saving to existing deterministic methods. The fastest existing methods generally [7] scale as \(O(M n^{4} \log n)\), or as \(O(M n^{3} \log n)\) in special cases [6].

The rest of the paper is organized as follows. In section 2, we introduce the Boltzmann equation and the LP scheme. In section 3, we review the TT format and the cross approximation algorithms. In section 4, we present the low-rank tensor method to compute the Boltzmann collision operator. In section 5, numerical results are presented to demonstrate the performance of the proposed method. Finally, we conclude the paper in section 6.

2 The space-homogeneous Boltzmann equation and its lifted version↩︎

In this section, we review the space-homogeneous Boltzmann equation and its properties. We then introduce a framework of evaluating the collision operator by a lifting-projection procedure, which was inspired by the theoretical breakthrough of Guillen and Silvestre [1] on well-posedness of the Landau equation. The same lifting technique was later adopted in the study of Boltzmann equation [2]. A detailed expository article on this topic in kinetic theory can be found in [3].

2.1 Boltzmann equation and its properties↩︎

The Boltzmann equation describes the behavior of a dilute gas of particles when the only interactions taken into account are binary elastic collisions. The space-homogeneous Boltzmann equation models the distribution function \(f(\mathbf{v},t)\), which represents the probability density of finding a particle with velocity \(\mathbf{v}\) at time \(t\). The general form of the 3D Boltzmann equation is given by \[\label{eq:boltzmann} \partial_{t}f(\mathbf{v},t) = q(f)(\mathbf{v},t)\tag{1}\] where \(q(f)\) is the collision operator, defined as \[q(f)(\mathbf{v},t) = \int_{\mathbb{R}^{3}}\int_{\mathbb{S}^{2}}\mathcal{B}(\mathbf{v}-\mathbf{w},\mathbf{\zeta})(f(\mathbf{v}')f(\mathbf{w}')-f(\mathbf{v})f(\mathbf{w}))\mathrm{d}\zeta\mathrm{d}\mathbf{w}.\] In this formula, \(\mathbf{v}\) and \(\mathbf{w}\) are the pre-collision velocities of two particles, \(\mathbf{v}'\) and \(\mathbf{w}'\) are the post-collision velocities, and \(\mathcal{B}(\mathbf{v}-\mathbf{w},\mathbf{\zeta})\) is the collision kernel, which depends on the relative velocity \(\mathbf{v}-\mathbf{w}\) and the scattering direction \(\mathbf{\zeta} \in \mathbb{S}^{2}\) (the unit sphere). These quantities are related by the following formulas: \[\begin{cases} \mathbf{v}' =& \frac{1}{2}\left(\mathbf{v} + \mathbf{w} + |\mathbf{v} - \mathbf{w}| \mathbf{\zeta}\right)\\ \mathbf{w}' =& \frac{1}{2}\left(\mathbf{v} + \mathbf{w} - |\mathbf{v} - \mathbf{w}| \mathbf{\zeta}\right). \end{cases}\]

In this work, for simplicity, we will focus on the variable hard sphere (VHS) model, where the collision kernel is given by \[\mathcal{B}(\mathbf{v}-\mathbf{w},\mathbf{\zeta}) = \frac{1}{4\pi}|\mathbf{v}-\mathbf{w}|^{\gamma},\] with \(\gamma\) being a parameter that characterizes the type of interaction between particles. For instance, \(\gamma = 1\) corresponds to hard spheres, while \(\gamma = 0\) corresponds to Maxwell molecules. Nevertheless, the numerical method we propose in this work can be extended to other collision kernels.

The Boltzmann equation has several important properties. Since the collision operator satisfies the following properties: \[\int_{\mathbb{R}^{3}}q(f)\begin{bmatrix} 1 \\ \mathbf{v}\\ |\mathbf{v}|^{2} \end{bmatrix}d\mathbf{v} = 0,\] i.e., the equation must conserve mass, momentum and energy: for any \(t > 0\), we have \[\int_{\mathbb{R}^{3}}f(\mathbf{v},t)\phi(\mathbf{v})d\mathbf{v} = \int_{\mathbb{R}^{3}}f(\mathbf{v},0)\phi(\mathbf{v})d\mathbf{v}, \forall~ \phi(\mathbf{v}) \in \{1, \mathbf{v}, |\mathbf{v}|^{2} \}.\] Moreover, testing Equation 1 with \(\log f\) and integrate over \(\mathbb{R}^{3}\), we have \[\begin{align} \frac{d}{dt} \int_{\mathbb{R}^{3}}f\log f d\mathbf{v} =& \int_{\mathbb{R}^{3}}q(f)\log f d\mathbf{v}\\ =& \int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\int_{\mathbb{S}^{2}}\mathcal{B}(\mathbf{v}-\mathbf{w},\mathbf{\zeta})(f(\mathbf{v}')f(\mathbf{w}')-f(\mathbf{v})f(\mathbf{w}))\log f(\mathbf{v})d\zeta d\mathbf{w}d\mathbf{v}\\ =& \frac{1}{4}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\int_{\mathbb{S}^{2}}\mathcal{B}(\mathbf{v}-\mathbf{w},\mathbf{\zeta})(f(\mathbf{v}')f(\mathbf{w}')-f(\mathbf{v})f(\mathbf{w}))\log \frac{f(\mathbf{v})f(\mathbf{w})}{f(\mathbf{v}')f(\mathbf{w}')}d\zeta d\mathbf{w}d\mathbf{v}\\ \leq& 0. \end{align}\] Hence the entropy defined as \(H(f) \mathrel{\vcenter{:}}= \int_{\mathbb{R}^{3}}f\log f d\mathbf{v}\) is non-increasing in time.

2.2 The scheme based on lifting-projection procedure↩︎

In [1], the authors resolved a long-standing open problem: the existence of global smooth solutions to the space-homogeneous Landau equation with Coulomb interaction. The key of the proof is to show decay of Fisher information, which was achieved through a lifting-projection procedure. The same lifting-projection technique was later extended to Boltzmann equation with very soft potentials [2]. Despite its huge impact on analysis, we have not seen existing work that adopts this technique in numerical methods. This paper is an effort in this direction.

To illustrate the idea of the lifting-projection procedure, we first introduce the concept of tangent flow as a generalization of the forward Euler scheme.

For a function \(F = F(\mathbf{v},\mathbf{w})\) defined on \(\mathbb{R}^{3} \times \mathbb{R}^{3}\), we define the lifted collision operator \(Q\) by \[Q(F)(\mathbf{v},\mathbf{w}) \mathrel{\vcenter{:}}= \int_{\mathbb{S}^{2}} \mathcal{B}(\mathbf{v}-\mathbf{w},\mathbf{\zeta}) \left(F(\mathbf{v}',\mathbf{w}') - F(\mathbf{v},\mathbf{w})\right)d\mathbf{\zeta}.\]

Let the projection operator \(\Pi\) be defined as \[\Pi F(\mathbf{v}) \mathrel{\vcenter{:}}= \int_{\mathbb{R}^{3}} F(\mathbf{v},\mathbf{w}) d\mathbf{w},\] and denote the tensor product of two functions as \[(f\otimes g) (\mathbf{v}, \mathbf{w})\mathrel{\vcenter{:}}= f(\mathbf{v})g(\mathbf{w}),\] then the Boltzmann collision operator \(q(f)\) satisfies: \[q(f) = \Pi Q (f\otimes f).\] On the time interval \((t_{0}, t_{0}+T)\), suppose that \(F(\mathbf{v},\mathbf{w},t)\) solve the lifted equation (a linear Kac master equation) \[\label{eq:2-kac} \partial_{t} F = Q (F)\tag{2}\] with initial condition \(F(\mathbf{v},\mathbf{w},t_{0})=f(\mathbf{v},t_{0})f(\mathbf{w},t_{0})\), then we have \[\left.\partial_{t} \Pi F(\mathbf{v},t) \right|_{t=t_{0}}= \Pi Q(f_{0}\otimes f_{0}) = \int_{\mathbb{R}^{3}} Q(f_{0}\otimes f_{0})d\mathbf{w} = q(f_{0}),\] which means the projection of the lifted solution \((t, \Pi F (\cdot, t))\) is a tangent flow of collision operator \(q\) at \((t_{0}, f(\cdot, t_{0}))\). We call it the lifting-projection(LP) flow associated with \(q\).

Given a discrete solution \(f_{h}^{(n)}\), we propose the following discrete LP scheme to obtain \(f^{(n+1)}_{h}\).

1. Solve the lifted Boltzmann equation in one time step \[\label{eq:discrete95lifted} F^{(s+1)}_{h} = f^{(s)}_{h}\otimes f^{(s)}_{h} + \sqint_{s\Delta t}^{(s+1)\Delta t } Q_{h}F_{h}(\tau) d\tau,\tag{3}\] where \(Q_{h}\) is the discretized version of lifted Boltzmann operator \(Q\), \(F_{h}(\tau)\) solves the equation \(\partial_{t} F_{h} = Q_{h} F_{h}\) on \((s \Delta t, (s+1)\Delta t)\) with initial condition \(f^{(s)}_{h}\otimes f^{(s)}_{h}\), and \(\sqint\) represents the numerical integration in time. In this paper, we only consider the forward Euler time integrator, i.e. \(\sqint_{s\Delta t}^{(s+1)\Delta t } Q_{h}F_{h}(\tau) d\tau = \Delta t Q_{h}F_{h}(s \Delta t)\)

2. Apply projection \[f^{(s+1)}_{h} =\Pi F^{(s+1)}_{h}.\]

2.3 Time integration of the Boltzmann equation↩︎

As have been shown in [22], the LP scheme introduced above is no more than first order accurate in time. To achieve higher order accuracy, we can adopt for example a high-order Runge-Kutta scheme. In this paper, we use the 4th order Runge-Kutta as follows: \[\label{eq:RK4} \begin{align} q_1 &= q\left(f_h^{(s)}\right),\\ q_2 &= q\left(f_h^{(s)} + \tfrac{\Delta t}{2}\,q_1\right),\\ q_3 &= q\left(f_h^{(s)} + \tfrac{\Delta t}{2}\,q_2\right),\\ q_4 &= q\left(f_h^{(s)} + \Delta t\,q_3\right),\\ f_h^{(s+1)} &= f_h^{(s)} + \frac{\Delta t}{6}\bigl(q_1 + 2q_2 + 2q_3 + q_4\bigr), \end{align}\tag{4}\] where \[\label{eq:subtract} q(g) = \Pi Q_h (g\otimes g)= \frac{1}{\Delta t} \left(\Pi (I+\Delta t Q_h)(g\otimes g) - g\right),\tag{5}\] which means we first solve the lifted equation 3 by forward Euler method, project back to 3D, and then subtract the original function to obtain the time differencing.

3 Review of the TT format↩︎

Although the lifted Boltzmann equation is a linear equation, its solution is a 6D function, which makes it infeasible to store and compute with the full tensor representation when the grid size is large [23]. To overcome the curse of dimensionality, we use the tensor train (TT) format [21] to represent the solution.

In this section, we review the TT format and the cross approximation algorithm which will be used to solve the lifted equation. The TT format is a low-rank representation of high-dimensional tensors, which allows us to efficiently store and compute with high-dimensional data. The cross approximation algorithms are used to compute the TT representation without ever computing or storing the full tensor. Instead, we only need to sample a small number of carefully chosen entries, which is essential for the efficiency of our low-rank Boltzmann solver.

3.1 Definition and operations↩︎

A tensor \(A \in \mathbb{R}^{n_{1} \times n_{2} \times \cdots \times n_{d}}\) is said to be in TT format [21] if its element can be represented as: \[\label{eq:TT95format} A(i_{1}, i_{2}, \cdots, i_{d}) = A^{(1)}_{\alpha_{0},\alpha_{1}}(i_{1})A^{(2)}_{\alpha_{1},\alpha_{2}}(i_{2})\cdots A^{(d)}_{\alpha_{d-1},\alpha_{d}}(i_{d})\tag{6}\] where \(A^{(\mu)} \in \mathbb{R}^{r_{\mu-1} \times n_{\mu} \times r_{\mu}}\) are the cores, and \(r_{\mu}\) are the ranks of the TT representation. When the ranks are low, the TT format offers tremendous saving for high-dimensional problems by reducing the storage cost from \(O(n^d)\) to \(O(dnr^2)\) where \(n\) is the maximal dimension of the tensor and \(r\) is the maximal rank of the cores. We briefly review some of the basic operations of the TT format, while refering the readers to [21] for more details.

Summation↩︎

The cores of the sum \(C = A + B\) are given by \[C^{(\mu)}(i_{\mu}) = \begin{bmatrix} A^{(\mu)}(i_{\mu}) & 0 \\ 0 & B^{(\mu)}(i_{\mu}) \end{bmatrix},~\mu = 2, \cdots, d-1,\] and \[C^{(1)}(i_{1}) = \begin{bmatrix} A^{(1)}(i_{1}) & B^{(1)}(i_{1}) \end{bmatrix},~C^{(d)}(i_{d}) = \begin{bmatrix} A^{(d)}(i_{d}) \\ B^{(d)}(i_{d}) \end{bmatrix}.\] As a result, the ranks of the sum are given by \(r_{\mu}^{C} = r_{\mu}^{A} + r_{\mu}^{B},~\mu = 2, \cdots, d-1\).

Contraction↩︎

The advantage of the TT format is that many operations can be implemented in a complexity proportional to \(n\) instead of \(n^{d}\). For example, consider the following one-dimensional contraction: \[C(i_{1}, \cdots, i_{\mu-1}, i_{\mu+1}, \cdots, i_{d}) = \sum\limits_{i_{\mu}=1}^{n_{\mu}} A(i_{1}, i_{2}, \cdots, i_{d})u(i_{\mu}).\] If \(A\) is in TT format as in Equation 6 , then the cores of \(C\) are given by \[C^{(1)} = A^{(1)}, \cdots, C^{(\mu-1)} = A^{(\mu-1)} \left(\sum\limits_{i_{\mu}=1}^{n_{\mu}} A^{(\mu)}(i_{\mu})u(i_{\mu}) \right)\in \mathbb{R}^{r_{\mu-2} \times r_{\mu}},~C^{(j-1)} = A^{(j)},~j = \mu+1, \cdots, d.\] Multidimensional contraction can be implemented in a similar way with complexity \(\mathcal{O}(nr^{2})\).

Rounding↩︎

After performing operations on tensors in TT format (e.g. summation), the TT ranks may increase. To control the ranks and maintain the efficiency of the representation, one needs to perform a rounding procedure, which can be seen as a generalization of truncated singular value decomposition (SVD). The rounding procedure allows us to approximate the tensor with a lower rank representation while controlling the approximation error. The rounding procedure [21] is listed in Algorithm 1. It was shown in [21] that the rounding procedure can be performed with complexity \(\mathcal{O}(dnr^{3}).\)

3.2 Cross approximation for TT↩︎

Cross approximation [24] is a class of algorithms for constructing a TT representation of a high-dimensional tensor accessing a few entries of the tensor, analogous to skeleton (CUR) decomposition for matrices. The key challenge is to select a small set of multi-indices at which to evaluate the tensor so that the resulting low-rank approximation is nearly optimal.

There are several cross approximation algorithms available in the literature [25], such as TT-Cross [24], the DMRG greedy cross approximation algorithm [26] etc. These algorithms differ in the way they select the indices for evaluation and update the cores, but they all aim to efficiently compute the TT representation with controlled ranks and approximation error. This is an active research area in multi-linear algebra. In this paper, we adopt the DMRG greedy cross method [26], which will be reviewed below.

The TT-cross approximation is given by: \[\label{eq:TTCA} \begin{align} &A(i_1,\dots,i_d)\\ \approx&\sum_{s_{1} \in \mathcal{I}^{\le 1},t_{1} \in \mathcal{I}^{>1}, \cdots,s_{d-1} \in \mathcal{I}^{\le d-1},t_{d-1} \in \mathcal{I}^{>d-1}} A(i_1,t_1) \left[ A(\mathcal{I}^{\le 1},\mathcal{I}^{>1}) \right]^{-1}(s_{1}, t_{1})\\ &A(s_1, i_2, t_{2})\left[ A(\mathcal{I}^{\le 2},\mathcal{I}^{>2}) \right]^{-1}(s_{2}, t_{2}) \cdots\cdots A(s_{d-1}, i_d) \\ =&\sum_{s_{k} \in \mathcal{I}^{\le k},t_{k} \in \mathcal{I}^{>k}, k = 0, 1, \cdots, d} \left(\prod_{k=1}^{d} A(s_{k-1}, i_k, t_{k}) \left[ A(\mathcal{I}^{\le k}, \mathcal{I}^{>k}) \right]^{-1}(s_{k}, t_{k})\right). \end{align}\tag{7}\] Here \(\mathcal{I}^{\le k}\) and \(\mathcal{I}^{>k}\) denote the positions of \(r_k\) rows and columns in the \(k\)-th unfolding \(A^{\{k\}}\mathrel{\vcenter{:}}= A(\overline{i_{1}\cdots i_{k}}, \overline{i_{k+1}\cdots i_{d}})\). To unify the notation, we introduce the border sets \(\mathcal{I}^{\le 0} = \{1\}\) and \(\mathcal{I}^{> d} = \{1\}\). The key step of the cross approximation algorithms is to select the indices \(\mathcal{I}^{\le k}\) and \(\mathcal{I}^{>k}\) for \(k = 1, 2, \cdots, d-1\). The selection of these indices determines the accuracy and efficiency of the approximation.

In [26], the authors propose a greedy cross approximation algorithm that selects the indices based on the maximum-volume principle, which aims to maximize the determinant of the submatrix formed by the selected indices. The algorithm considers interpolation sets \(\mathcal{I}^{\le k}\) and \(\mathcal{I}^{>k}\) that satisfy the following nestedness condition: \[\label{eq:nestedness} \begin{align} \overline{i_{1}\cdots i_{k}} \in \mathcal{I}^{\le k} &\Rightarrow \overline{i_{1}\cdots i_{k-1}} \in \mathcal{I}^{\le k-1}\\ \overline{i_{k+1}\cdots i_{d}} \in \mathcal{I}^{>k} &\Rightarrow \overline{i_{k+2}\cdots i_{d}} \in \mathcal{I}^{>k+1} \end{align}\tag{8}\] We describe the algorithm in Algorithm 2 for completeness. The stopping criterion is that the relative \(l_2\) error on a random test set is less than the given tolerance or the number of sweep reaches given threshold. This algorithm requires \(\mathcal{O}(dnr^{2}X)\) evaluation of tensor elements and \(\mathcal{O}(dnr^{3}X)\) additional operations, with \(X\) being the number of sweeps.

4 Low-rank Boltzmann solver based on TT↩︎

With the TT format and the cross approximation algorithms, we are now ready to describe the fast scheme for the Boltzmann equation based on TT. For simplicity, we only describe the scheme for the 3D case. The 2D case is simpler, and the description is skipped for brevity.

4.1 Discretization↩︎

As a standard practice, we need to choose a cut-off for the velocity domain. For the 6D veclocity space, we choose \((\mathbf{v},\mathbf{w}) \in (-L, L)^{3} \times (-L, L)^{3}\) as the truncated domain. As have been illustrated in Remark 3, on each sphere \(\mathcal{B}(\frac{z}{2},\lvert \mathbf{u} \rvert)\), the lifted Boltzmann drives the lifted solution toward spherical average through a reaction equation \(F_{t} = \lvert \mathbf{u} \rvert^{\gamma} \left(\mathcal{A}F - F\right)\). When calculating the spherical average, we assume that the lifted solution \(F\) is equal to zero outside of truncated domain.

We discretize the 3D truncated domain \((-L, L)^3\subset \mathbb{R}^{3}\) with a uniform grid of size \(n\) in each dimension, and denote the grid points as \(\{(v_{x}(i_{1}), v_{y}(i_{2}), v_{z}(i_{3}))\}\), where \(i_{1}, i_{2}, i_{3} = 1, 2, \cdots, n\). In this paper, we assume the number of points in different directions to be the same, though the method is easily applicable with different number of points in different directions. The grid points are given by \[\label{eq:gridpoint} \left(v_{x}(i_{1}), v_{y}(i_{2}), v_{z}(i_{3})\right) = \left(-L + (i_{1}-\frac{1}{2})h, -L + (i_{2}-\frac{1}{2})h, -L + (i_{3}-\frac{1}{2})h\right).\tag{9}\] where \(h = \frac{2L}{n}\) is the grid size. Note that for simplicity we assume the numbers of nodes are the same in all directions. The discrete solution is denoted as \(f_{h}(\mathbf{v})\), which is a tensor of size \(n \times n \times n\): \[A(i_{1}, i_{2}, i_{3}) = f_{h}\left(v_{x}(i_{1}), v_{y}(i_{2}), v_{z}(i_{3})\right).\] The lifted solution is defined on the Cartesian product of the velocity grid with itself, and its discrete version is denoted as \(F_{h}(\mathbf{v}, \mathbf{w})\), which is a tensor of size \(n \times n \times n \times n \times n \times n\): \[\label{eq:xyordering} G(i_{1}, i_{2}, i_{3}, i_{4}, i_{5}, i_{6}) = F_{h}\left(v_{x}(i_{1}), v_{y}(i_{3}), v_{z}(i_{5}), w_{x}(i_{2}), w_{y}(i_{4}), w_{z}(i_{6})\right).\tag{10}\] We call this ordering the \(xyz\)-ordering, which is the recommended ordering from our numerical experiment since it fully exploits the low-rank property between the velocity dimensions.

When solving the VHS model, the key step is the computation of spherical averaging operator \[\mathcal{A} F(\lvert \mathbf{u} \rvert, \sigma;\mathbf{z}) \mathrel{\vcenter{:}}= \frac{1}{4\pi}\int_{\mathbb{S}^{2}} F(\lvert \mathbf{u} \rvert, \zeta;\mathbf{z}) d\zeta,\] where \(\mathbf{u} = \mathbf{v} - \mathbf{w}\), \(\mathbf{z} = \mathbf{v} + \mathbf{w}\), and \(\mathbf{\sigma} = \mathbf{u}/|\mathbf{u}|\).

Denote the spherical quadrature points and weights as \(\left\{(\mathbf{\zeta}_{i}, q_{i})\right\} \subset \mathbb{S}^{2} \times \mathbb{R}\), where \(i = 1, 2, \cdots, M\), we define \[A_{h} F_{h}(\mathbf{v},\mathbf{w}) = \sum\limits_{i=1}^{M} q_{i} F_{h}(\frac{\mathbf{z} + |\mathbf{u}|\mathbf{\zeta}_{i}}{2}, \frac{\mathbf{z} - |\mathbf{u}|\mathbf{\zeta}_{i}}{2}).\] Tensor product of one dimensional quadrature rules along the azimuthal and polar directions can be computationally expensive because the number of quadrature points grows quadratically with respect to the number of points in each direction. To address this issue, we use the Lebedev quadrature rule [27], which is a family of quadrature rules for integration on the sphere that achieves high accuracy with a relatively small number of quadrature points.

Since quadrature points on the sphere does not necessarily lie on the discrete velocity grid, special treatment is needed. When it falls outside of the truncated domain \((-L, L)^{3}\), we simply use zero padding. For points inside the box \((-L,L)^{3}\), we use one of the following interpolations:

• A 4-point stencil for interpolation in each dimension, which has order 4 accuracy. The interpolation weights are computed using the Lagrange polynomial basis.

• Spectral interpolation with a global stencil using FFT.

More details will be described below.

4.2 First-order in time scheme↩︎

Suppose we have the solution \(f^{(s)}_{h}\) in TT format. To obtain the updated solution \(f^{(s+1)}_{h}\) with first order accuracy in time, we take the following steps:

1. Compute the lifted solution \(F_{h}^{(s+1)}\) in TT format using cross approximation (Algorithm 6).

2. Project the lifted solution \(F_{h}^{(s+1)}\) to the original velocity space to obtain \(f_{h}^{(s+1)}\) in TT format.

3. Perform TT-rounding on \(f_{h}^{(s+1)}\) to control the ranks of the solution.

These steps are summarized in Algorithm 5, where a simple flowchart is plotted in Figure 3.

Specifically, the cross approximation of the 6D lifted tensor is done by Greedy DMRG (Algorithm 2). The stopping criterion for the Greedy DMRG is that the relative \(l_2\) error on a random test set is less than the given tolerance or the number of sweep reaches given threshold. In practice, we set the max number of sweep to be the same as grid size in each dimension. Moreover, the tolerance parameter is set to be on the order of the local truncation error, for example, \(\frac{1}{2} \Delta t (h^{4} + (\Delta t)^{4})\) for cubic interpolation with RK4 scheme. Point evaluations for the cross approximations are described in Algorithm 6, which is plotted in Figure 4 (for the 2D case with cubic interpolation for simplicity). The details of cubic and spectral interpolations are described in Algorithms 7 and 8. For cubic interpolation (Algorithm 7), the computational cost scales as \(\mathcal{O}(Mr^{2})\). For spectral interpolation (Algorithm 8), the computational cost scales as \(\mathcal{O}(Mr^{2}n)\). Finally, the purpose of the TT-rounding is to limit the rank of the TT. We suggest using the local truncation error as the rounding tolerance, for example, \(\frac{1}{20} \Delta t (h^{4} + (\Delta t)^{4})\) for cubic interpolation with RK4 scheme.

4.3 Higher-order in time and conservation correction↩︎

The extensions to higher-order in time, e.g. RK4 is described in Algorithm 9. The basic idea is to embed the first-order algorithm with a Runge-Kutta scheme. Here the tolerance of the TT-rounding is set to be on the order of the local truncation error of the scheme. In practice we are using \(\frac{1}{20} \Delta t (h^{4} + (\Delta t)^{4})\) for cubic interpolation with RK4 scheme. For numerical tests with spectral interpolation, we set the TT-rounding and cross approximation tolerances to a small enough constant that can never be reached, which is \(10^{-10}\) in our numerical experiments.

In addition, we develop a conservation correction step that is friendly to the TT format. Due to numerical discretization, domain truncation and cross approximation, the conservation of mass, momentum and energy is generally lost. We propose to use a Lagragian multiplier approach and solve the following constrained optimization problem to compute the correction term \(\mu g_{h}\) to be added to the computed solution \(f_{h}\): \[\min_{g_{h}} \lVert g_{h}\rVert_{L^{2}_{\mu}}^{2} \text{ s.t. } (g_{h}, \phi_{k})_{L^{2}_{\mu}} = \mathcal{M}_{k} - (f_{h}, \phi_{k})_{L^{2}}, \quad k = 1, 2, \cdots, 5\] where \(\{\phi_{k}\} = \{1, v_{x}, v_{y}, v_{z}, |\mathbf{v}|^{2}\}\) to enforce mass, momentum and energy conservation [8]. The weight \(\mu\) is set to be \(\mu = \exp (-\lvert \mathbf{v}\rvert^{2}/2\sigma^{2})\), in practice we choose \(\sigma = 1\). (When \(\mu \equiv 1\) one recovers the correction in [8].) We choose this weight to accommodate the zero boundary conditions on the boundary, and in practice, it shows better performance in accuracy.

Let \(A^{n^{3} \times 1}\) be the vectorization of the computed solution \(f_{h}\), \(M^{5 \times 1}\) be the conserved quantity, \(W^{n^{3}\times n^{3}}\) be the matrix associated to the weight \(\mu\), and \(\Phi^{n^{3}\times 5}\) be the matrix whose columns are the vectorization of \(\{\phi_{k}\}\), then the optimization problem can be rewritten as \[\min_{G} G^{T}WG \text{ s.t. } \Phi^{T} (W G + A) = M.\] Analogous to [8], the constrained optimization problem can be solved explicitly as follows: \[G = \Phi(\Phi^{T}W\Phi )^{-1}( M - \Phi^{T}A ).\]

In implementation, this correction is performed in the TT format. Since \(\{\phi_{k}\}\) are low-rank functions in the velocity space, \(\Phi\) can be represented in TT format with low rank. Consequently, the correction term \(WG\) can be computed with complexity \(\mathcal{O}(dnr^{2})\) instead of \(\mathcal{O}(n^{3})\) for the full tensor. In addition, we note the correction is of low-rank in TT format. This is because \(g_{h}\) is always a linear combination of \(\{\phi_{k}\} = \{1, v_{x}, v_{y}, v_{z}, \mathbf{v}^{2}\}\). If we define \[\begin{align} &g_{1}(v_{x},1) = 1,\;g_{1}(v_{x},2) = \kappa_{1}v_{x} + \kappa_{4}v_{x}^{2},\\ &g_{2}(v_{y},1,1) = 1 + \kappa_{2}v_{y} + \kappa_{4}v_{y}^{2},\;g_{2}(v_{z},2,2) = 0, g_{2}(v_{y}, 1,2) = g_{2}(v_{y}, 2,1)=1,\\ &g_{3}(v_{z},1) = 1,\;g_{3}(v_{z},2) = \kappa_{3}v_{z} + \kappa_{4}v_{z}^{2}, \end{align}\] then \[1 + \kappa_{1}v_{x} + \kappa_{2}v_{y} + \kappa_{3}v_{z} + \kappa_{4} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) = \sum_{\alpha_{1} = 1}^{2}\sum_{\alpha_{2} = 1}^{2} g_{1}(v_{x},\alpha_{1})g_{2}(v_{y},\alpha_{1}, \alpha_{2})g_{3}(v_{z},\alpha_{2}).\] and it follows that \[\begin{align} &\left[1 + \kappa_{1}v_{x} + \kappa_{2}v_{y} + \kappa_{3}v_{z} + \kappa_{4} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2})\right]\mu_{1}(v_{x})\mu_{2}(v_{y})\mu_{3}(v_{z}) \\ =& \sum_{\alpha_{1} = 1}^{2}\sum_{\alpha_{2}=1}^{2} \mu_{1}g_{1}(v_{x},\alpha_{1})\mu_{2}g_{2}(v_{y},\alpha_{1}, \alpha_{2})\mu_{3}g_{3}(v_{z},\alpha_{2}). \end{align}\] Therefore, the correction term \(\mu g_{h}\) is always a TT with max rank equal to 2.

4.4 Complexity analysis↩︎

Suppose that the lifted 6D solution has TT-ranks \(\{R_{1}, \cdots, R_{5}\}\) with maximum \(R\), and the projected 3D solution has TT-ranks \(\{r_{1},r_{2}\}\) with maximum \(r\).

The cost of TT cross approximation is dominant in Algorithm 5. The greedy cross approximation Algorithm 2 takes \(\mathcal{O}(nR^{2}X)\) evaluations of the tensor elements, and \(\mathcal{O}(nR^{3}X)\) flops for other operations, with \(X\) being the number of sweeps. It takes \(\mathcal{O}(Mr^{2})\) interpolations for each query of Algorithm 6. Each cubic interpolation (Algorithm 7) takes \(\mathcal{O}(1)\) flops and each spectral interpolation (Algorithm 8) takes \(\mathcal{O}(n)\) flops. Therefore the total cost of one time step is in the order of \(\mathcal{O}(nMmR^{2}r^{2}X + nR^{3}X)\), where \(m\) is the width of interpolation stencil.

Finally, we would like to briefly comment on alternative implementations based on TT. We have conducted numerical investigations into two alternative approaches, and below is our findings.

• (Alternative 1) Perform cross approximation directly on the lifted operator \(Q(g\otimes g)\) rather than \(\left(I+ \Delta t Q\right)(g\otimes g)\). The solver does not perform well. This is because when \(g\) is almost in equilibrium, \(Q(g \otimes g )\) is near zero, which is not computational robust for cross approximation.

• (Alternative 2) Directly compute the original Boltzmann equation 1 with TT without lifting-projection. If we perform a TT cross approximation on the forward Euler scheme for 1 , this will not take advantage of fast operations of the lifted solution. For example, the cross approximation Algorithm 2 needs \(\mathcal{O}(nr^{2})\) point evaluations of tensor elements.

  For each 3D velocity \(\mathbf{v}\), evaluation of \(\Pi(I + Q \Delta t) (g\otimes g)(\mathbf{v})\) requires \(\mathcal{O}(M n^{3})\) quadrature points since \(\Pi Q\) includes integration on \(\mathbb{R}^{3}\times \mathbb{S}^{2}\). If \(g\) is given as a 3D TT with max rank \(r\), then the total cost will be \(\mathcal{O}(Mn^{4}r^{4})\), considering that the cubic interpolation of \(g\) at off-grid velocity takes \(\mathcal{O}(r^{2})\) flops. In most cases this is much more expensive than the approach we propose.

5 Numerical results↩︎

In this section, we conduct numerical simulations on 2D and 3D benchmark problems, focusing on the performance of the method in terms of accuracy, speed and low-rank properties. In what follows, we choose \(L=6.4,\) i.e. using \((-6.4, 6.4)^{2}\) or \((-6.4, 6.4)^{3}\) as the truncated domain.

5.1 Experiments on 2D Boltzmann equation↩︎

We use the 2D case to study the effect of the ordering for the ranks and behaviors of the solution, as well as benchmarking on numerical accuracy.

5.1.1 Ordering of variables for TT↩︎

As has been shown in Section 4, the complexity of the low-rank method depends on the rank of TTs. Therefore, it is crucial to investigate how the ordering of variables affects the TT rank of the lifted solution. We consider:

1. The \(vw\)-ordering \((v_{x}, v_{y}, w_{x}, w_{y})\) (as defined in ?? ).

2. The \(xy\)-ordering \((v_{x}, w_{x}, v_{y}, w_{y})\) (as defined in 10 for the 3D case).

The orderings are tested on three different initial conditions with the Maxwell molecules where \(\gamma = 0\):

1. The BKW [28], [29] reference solution: \[\label{eq:bkw2d} f_{0}(\mathbf{v}) = |\mathbf{v}|^{2}\exp(-|\mathbf{v}|^{2})\tag{11}\]

2. Two Gaussians centered at \(\mathbf{v}_{1} = (-2, 0)\) and \(\mathbf{v}_{2} = (2, 0)\): \[\label{eq:2gaussian2d} f_{0}(\mathbf{v})=\exp(-\frac{|\mathbf{v}-\mathbf{v}_{1}|^2}{2}) + \exp(-\frac{|\mathbf{v}-\mathbf{v}_{2}|^2}{2})\tag{12}\]

3. Anisotropic Gaussian centered at the origin: \[\label{eq:anigaussian2d} f_{0}(\mathbf{v}) = \exp(-\frac{(v_{x} + v_{y})^{2}}{8}-\frac{(v_{x} - v_{y})^{2}}{2})\tag{13}\]

We note that the last example is more challenging for low-rank methods due to its intrinsic higher rank between \(v_{x}\) and \(v_{y}\).

In each test, let the mesh size \(h\) be \(\frac{2L}{n} = \frac{2\times6.4}{32} = 0.4\). For purpose of benchamarking, we use full tensor, i.e. in Algorithm 5, we keep the same interpolation and quadrature rules, but replace the TT cross-approximation step with full-tensor construction. We compute the full 4D tensor associated to the lifted solution at the first time step \(t = \Delta t = 0.1\), denoted by \(G\), and perform TT-SVD to obtain the optimal low-rank approximation \(G_{*}\) such that: \[\left\lVert G - G_{*} \right\rVert_{F} \leq \varepsilon \left\lVert G \right\rVert_{F},\] where \(\varepsilon\) is set to be \(10^{-5}\).

The results are summarized in Table 1, where we report the TT ranks of the lifted solution \(G\) with respect to various orderings and initial conditions. It is clear the advantage of the \(xy\)-ordering compared to \(vw\)-ordering for BKW and two Gaussian test cases where the true solution exhibits low-rank properties. For the anisotropic Gaussian case, the two methods of ordering show similar ranks. Since the main use case for low-rank methods are when the solution \(f\) exhibits low-rank properties, the \(xy\)-(or \(xyz\)-ordering in 3D) is well motivated.

5.1.2 Accuracy test↩︎

To verify the order of accuracy, we present the results of the following test: Using maxwell molecule model, i.e. \(\mathcal{B}(\mathbf{v}-\mathbf{w},\zeta) = \frac{1}{2\pi}\), we let the BKW initial condition \[f(\mathbf{v}, t) = \frac{1}{2\pi K(t)^{2}}\exp\left(-\frac{|\mathbf{v}|^{2}}{2K(t)}\right)\left( 2K(t)-1 + \frac{1-K(t)}{2K(t)}|\mathbf{v}|^{2}\right),\] where \(K(t) = 1 - 0.5\exp{(-t/8)}\), evolve from \(t=0\) to \(t=4\), and compute the relative \(L^{2}\) error by comparing the computed solution with the reference solution given by the BKW formula. We use RK4 as time integrator with cubic interpolation. The run uses \(M=64\) circular quadrature points, cross-approximation tolerance \(\frac{1}{2}\Delta t(h^4+(\Delta t)^4)\), at most \(n\) greedy sweeps, and TT-rounding tolerance \(10^{-16}\).

The results are shown in Table 2. We can clearly observe 4th order convergence in space and time.

5.2 Experiments on 3D Boltzmann equation↩︎

5.2.1 Accuracy test on cubic interpolation↩︎

First we test the RK4 method: Algorithm 9 with cubic interpolation: Algorithm 7. In these tests we use \(M=74\) Lebedev points, cross-approximation tolerance \(\frac{1}{2}\Delta t(h^{4}+(\Delta t)^4)\), a maximum number of greedy sweeps equal to the number of grid points in each velocity direction, and TT-rounding tolerance equal to one tenth of the cross-approximation tolerance.

We expect the overall order of accuracy of the proposed method to be 4. To verify the order of accuracy, we present the results of two experiments:

• Maxwell molecule model, i.e. \(\mathcal{B}(\mathbf{v}-\mathbf{w},\zeta) = \frac{1}{4\pi}\), we consider the BKW solution \[\label{eq:bkw3d} f(\mathbf{v}, t) = \frac{1}{(2\pi K(t))^{3/2}}\exp\left(-\frac{|\mathbf{v}|^{2}}{2K(t)}\right)\left( \frac{5K(t)-3}{2K(t)} + \frac{1-K(t)}{2K(t)^{2}}|\mathbf{v}|^{2}\right),\tag{14}\] where \(K(t) = 1 - \exp{(-(t+5.5)/6)}\), evolve from \(t=0\) to \(t=2\), and compute the relative \(L^{2}\) error by comparing the computed solution with the reference solution given by the BKW formula.

• Hard sphere model, i.e. \(\mathcal{B}(\mathbf{v}-\mathbf{w},\zeta) = \frac{1}{4\pi} \left\vert\frac{\mathbf{v}-\mathbf{w}}{6.4}\right\vert\), we let the two Gaussians initial condition \[\label{eq:twoG3d} f_{0}(\mathbf{v})=\exp(-\frac{|\mathbf{v}-\mathbf{v}_{1}|^2}{2}) + \exp(-\frac{|\mathbf{v}-\mathbf{v}_{2}|^2}{2})\tag{15}\] where \(\mathbf{v}_{1} = (2, 0, 0)\) and \(\mathbf{v}_{2} = (0, -2, 0)\), evolve from \(t=0\) to \(t=2\). Since there is no exact reference solution, we only perform self-convergence test.

The results are shown in Tables ¿tbl:tab:accuracy95cubic95mm? and [tbl:tab:accuracy95cubic95hs]. In both cases, we can clearly observe 4th order convergence upon mesh refinement.

In addition, the initial and final states are plotted in Figures 10 (c) and 10, which benchmark well with known solutions in the literature.

5.2.2 Accuracy test on spectral interpolation↩︎

Next we test the RK4 method: Algorithm 9 with spectral interpolation: Algorithm 8.

We expect the method to have spectral accuracy in velocity and 4th order accuracy in time. We solve the same test cases as in 14 and 15 . For the spectral runs we use \(M=74\) Lebedev points, cross-approximation tolerance \(10^{-10}\), at most \(n\) greedy sweeps, and TT-rounding tolerance \(10^{-10}\).

The results of Maxwell molecules with BKW initial condition are shown in Table ¿tbl:tab:accuracy95mm95spectral?. As we can see, extremely accurate numerical solution can be obtained with \(n=32,\) reflecting spectral accuracy. The error saturates at the level of \(10^{-8}\) due to the prescribed truncation tolerance and possibly the temporal discretization error.

The results of hard spheres with “two-Gaussians" initial condition are shown in Table [tbl:tab:accuracy95hs95spectral]. In this case, the solution is not strictly low-rank. We do not observe spectral accuracy because the cross approximation algorithm reaches max number of sweep and fails to meet the tolerance.

5.2.3 Complexity of cubic interpolation↩︎

Recall that with cubic interpolation, the complexity of the proposed method is \(\mathcal{O}(nMR^{2}r^{2}X + n R^3 X)\), where \(n\) is the number of grid points in each direction, \(M\) is the number of quadrature points on the sphere, \(R\) is the max TT rank of the 6d lifted solution \(F\), \(r\) is the TT-rank of the 3D projected solution \(f\), and \(X\) is the number of sweeps. In what follows we present the computational complexity with two test cases, 14 and 15 . The timing runs use \(M=74\) Lebedev points, a maximum number of greedy sweeps equal to the number of grid points per velocity direction, and TT-rounding tolerance equal to one tenth of the cross-approximation tolerance.

Figure 11 shows the maximal TT ranks of the numerical solutions over time for two test cases. We can see that the TT ranks of the lifted solutions \(R\) are higher than \(r\) and vary with initial conditions and time.

To show the performance of the proposed method, we let the system evolve from \(t=0\) to \(t=0.5\) with RK4 time discretization, where \(\Delta t = 0.1\). The average computational time per time step is shown in Tables 3 and 4 for the maxwell molecule model and hard sphere model respectively. The time cost is measured on MacBook Pro with M4 chip, with a MATLAB implementation. In particular, for tensor train operations we use the implementation provided in the GitHub repository by Oseledets et. al. [30], with slight modifications.

For fixed tolerance, the time cost is proportional to \(n\), since the rank of the lifted solution is almost independent of \(n\). For fixed \(n\), the time cost increases with respect to the accuracy requirement, since higher accuracy requires higher rank for the lifted solution. The BKW solution with Maxwell molecule is a special case where the lifted solution is strictly low-rank, which explains the low computational cost in Table 3. The computational results verify the complexity scaling. In particular, we note the importance of choosing appropriate tolerance according to local truncation error. If the tolerance is set to be too small, the cross approximation routine will fail to meet the tolerance and reach max number of sweep, which explains high computational cost.

5.2.4 Complexity of spectral interpolation↩︎

Recall that with spectral interpolation, the complexity of the proposed method is \(\mathcal{O}(n^{2}MR^2 r^2 X + nR^3X)\), where \(n\) is the number of grid points in each direction, \(M\) is the number of quadrature points on the sphere, \(R\) is the max TT rank of the 6d lifted solution \(F\), \(r\) is the TT-rank of the 3D projected solution \(f\), and \(X\) is the number of sweeps.

Following the same procedure and using the same parameter as in Section 5.2.3, the average time cost is shown in Tables 5 and 6. For the Maxwell molecule BKW test, the measured times are nearly independent of the tolerance because the lifted solution remains very low rank. For the hard sphere test, tighter tolerances increase the ranks and sweep counts, so the time grows with both \(n\) and the requested accuracy, consistent with the rank-dependent complexity estimate. When comparing Table 6 with Table 4, it can be observed that for \(n = 32\), the spectral solver is slightly faster than the cubic one, while in the \(n=128\) case, it is much slower. This observation is also consistent with the complexity analysis: when \(n\) is small, the auxiliary operation cost \(\mathcal{O}(nR^3X)\) in TT cross approximation dominates, and as \(n\) increases, the interpolation cost \(\mathcal{O}(n^{2}MR^2 r^2 X)\) becomes the dominant term since it scales quadratically with \(n\).

5.3 Necessity of conservation correction↩︎

In prior subsections, all the low-rank numerical results are obtained with conservation correction. In this subsection, we run the numerical tests with and without conservation correction to validate its effect. Using cubic interpolation and RK4, we set \(L=6.4\), \(M=74\) Lebedev points, \(n=64\), \(\Delta t=0.2\), a maximum number of greedy sweeps equal to the number of grid points per velocity direction, and TT-rounding tolerance equal to one tenth of the cross-approximation tolerance.

The numerical results are gathered in Figures 12 and 13. It is clear that the numerical results without conservation correction present large errors in the moments. Indeed, an essential condition for the lifting-projection procedure to work is that the 3D solution \(f\) should be a probability distribution s.t. \(\int_{\mathbb{R}^{3}} f(\mathbf{v}) d \mathbf{v} = 1\). More specifically, in Algorithm 9, when computing \[q(f^{(s)})= \Pi Q (f_{h}^{(s)}\otimes f_{h}^{(s)}) = \frac{\Pi (I + Q \Delta t) (f_{h}^{(s)}\otimes f_{h}^{(s)}) - f_{h}^{(s)}}{\Delta t},\] we have implicitly assumed that \[\Pi (f_{h}^{(s)}\otimes f_{h}^{(s)}) - f_{h}^{(s)} =0,\] which is only true when mass conservation is preserved in every step. The computational results validate the essential roles of conservation corrections.

6 Conclusion↩︎

We presented a novel conservative, low-rank TT Boltzmann solver consisting of two independent contributions. The first is the lifting-projection (LP) scheme, which is a computational framework motivated by recent theoretical breakthrough [1]–[3] that lifts the nonlinear 3D Boltzmann equation to a 6D linear Kac master equation, evolve for one time step, and projects the result back to 3D. The second contribution is a fast low-rank tensor method. We represent the lifted solution in TT format and compute it via a TT cross algorithm. With cubic interpolation, the cost of the proposed method is \(\mathcal{O}(nMR^{2}r^{2}X + n R^3 X)\), where \(n\) is the number of grid points in each direction, \(M\) is the number of quadrature points on the sphere, \(R\) is the max TT rank of the 6d lifted solution \(F\), \(r\) is the TT-rank of the 3D projected solution \(f\), and \(X\) is the number of sweeps. With spectral interpolation, the cost is \(\mathcal{O}(n^{2}MR^{2}r^{2}X)\). When the solution exhibits low-rank property, this method, which scales linearly with \(n,\) offers huge computational savings compared to existing methods in the literature. We further propose a simple conservation correction procedure that is compatible with TT format that ensures mass, momentum and energy conservation.

We perform extensive numerical tests, benchmarking the performance of the methods in terms of accuracy and computational efficiency. The numerical experiments demonstrate that the \(xyz\)-ordering in TT leads to significantly lower ranks for the lifted solution compared to the \(vw\)-ordering. Convergence tests show that the methods achieve the expected accuracy and numerical conservations are maintained.

Several directions remain open for future work. The method currently handles only the spatially homogeneous Boltzmann equation. Extending it to the full spatially inhomogeneous case is a natural next step. We also plan to investigate applications to Landau equation and study of the asymptotic preserving schemes.

Acknowledgements↩︎

This material is based upon work supported by the National Science Foundation under Grant No. DMS-2424139 while the authors were in residence at the Simons Laufer Mathematical Sciences Institute in Berkeley, California, during the Fall 2025 semester.

References↩︎

---

1. Department of Mathematics, Virginia Tech, Blacksburg, VA 24061 U.S.A. (kunhuang@vt.edu)↩︎

2. Department of Mathematics, Virginia Tech, Blacksburg, VA 24061 U.S.A. (yingda@vt.edu). Research is supported by DOE grant DE-SC0023164, AFOSR grant FA9550-25-1-0154 and Virginia Tech.↩︎

3. Department of Mathematics and Oden Institute for Computational Sciences and Engineering, University of Texas at Austin (gamba@math.utexas.edu). Research is supported by NSF grant DMS 2408263.↩︎