embodies the kinetic contribution[@condiff_transport_1969] due to translational, rotational and coupled velocities. On the other hand, the heat flux only shows a non-equilibrium contribution and, given the symmetry of $f_{leq}^{(1)}$, $$\begin{aligned}
 \boldsymbol{\mathfrak{Q}}({\boldsymbol{\mathfrak{r}}};t)=\frac{1}{2}\int d{\boldsymbol{\mathfrak{p}}}\,|\mathbf{c}|^2\,{\boldsymbol{\mathfrak{m}}}^{-1}\mathbf{c}\,f^{(1)}({\boldsymbol{\mathfrak{r}}},{\boldsymbol{\mathfrak{p}}};t)
=\frac{1}{2}\int d{\boldsymbol{\mathfrak{p}}}\,|\mathbf{c}|^2\,{\boldsymbol{\mathfrak{m}}}^{-1}\mathbf{c}\,f_{neq}^{(1)}({\boldsymbol{\mathfrak{r}}},{\boldsymbol{\mathfrak{p}}};t)\label{eq:24}.
 \end{aligned}$$ Now, making use of the identity $$\begin{aligned}
{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}}\cdot\left[({\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t)\otimes{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t))\rho({\boldsymbol{\mathfrak{r}}};t)\right] =
{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t){\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}}\cdot\left({\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t)\rho({\boldsymbol{\mathfrak{r}}};t)\right)+\rho({\boldsymbol{\mathfrak{r}}};t)\left({\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t)\cdot{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}}\right){\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t),\label{eq:25}\end{aligned}$$ in Equations (\[eq:23\]) and (\[eq:19\]), we arrive at the time-evolution equation $$\begin{aligned}
 \rho({\boldsymbol{\mathfrak{r}}}_1;t)\mathcal{D}_t{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_1;t)\label{eq:26}+&{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\cdot {\overline{\overline{\boldsymbol{\Pi}}}}_{neq}({\boldsymbol{\mathfrak{r}}}_1;t)+{\boldsymbol{\mathfrak{m}}}^{-1}\rho({\boldsymbol{\mathfrak{r}}}_1;t){\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}_1}\frac{\delta \mathcal{F}[\rho]}{\delta \rho({\boldsymbol{\mathfrak{r}}}_1)}
 +({\boldsymbol{\Upsilon}}_{\parallel}+{\boldsymbol{\Upsilon}}_{\bot}{\boldsymbol{\epsilon}})\rho({\boldsymbol{\mathfrak{r}}}_1;t){\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}}_1;t)\\
 +&N{\boldsymbol{\mathfrak{m}}}^{-1}\sum_{j=1}^{N}\int d{\boldsymbol{\mathfrak{r}}}^{N-1}d{\boldsymbol{\mathfrak{p}}}^N\,\widetilde{{\boldsymbol{\Upsilon}}}_{1j}({\boldsymbol{\mathfrak{r}}}^N){\boldsymbol{\mathfrak{p}}}_j{f^{(N)}}({\boldsymbol{\mathfrak{r}}}^N,{\boldsymbol{\mathfrak{p}}}^N;t) = 0\notag\end{aligned}$$ where we define the material derivative as $\mathcal{D}_t=(\partial_t+{\boldsymbol{\mathfrak{v}}}({\boldsymbol{\mathfrak{r}}};t)\cdot{\boldsymbol{\nabla}}_{{\boldsymbol{\mathfrak{r}}}})$ and ${\boldsymbol{\epsilon}}=\text{adiag}(m^{-1}\mathbb{I},m\mathbb{I}^{-1})$. At this stage, there are two terms which still have to be treated, namely the non-equilibrium and the friction tensor terms. The former will be analysed in what follows and the latter in Section \[subsec:enskog-closure\].

To treat the non-equilibrium term, we follow Grad’s method[@grad_kinetic_1949; @liboff_kinetic_2003]. Specifically, the one-body distribution is rewritten as $f^{(1)}=f^{(1)}_{leq}\,\left[1+\varPsi
\right]$ and $\Psi({\boldsymbol{\mathfrak{r}}},{\boldsymbol{\mathfrak{p}}};t)$ is expanded in a basis of orthogonal polynomials of the $\mathbf{p}\oplus\mathbf{L}$ space. In particular, we choose a basis of the multiple generalised Hermite polynomials[@dunkl_classical_2001] on $\mathbb{R}^d$, $$\begin{aligned}
P^{[\alpha]}_{n,\boldsymbol{\ell}}(x_1,\dots,x_d)=H_{\ell_1}^{[\alpha]}(x_1)\dots H_{\ell_d}^{[\alpha]}(x_d)\label{eq:27},
\quad \sum\ell_k=n\end{aligned}$$ with $n\in\mathbb{N}_0=\mathbb{Z}^+\cup\{0\}$, $\boldsymbol{\ell}\in\mathbb{N}^d_0$ and $H_{k}^{[\alpha]}(x)\equiv
H_{k}(x/\sqrt{\alpha})/\sqrt{k!}$ the generalised Hermite polynomials, which are an orthogonal basis for $L^2\left(\mathbb{R}\right)$ with respect to the weight $(2\pi\alpha)^{-1/2} \exp( -x^2/(2\alpha) )$. Denoting $\mathbf{e}_j=(\delta_{jk})$, the first $P^{[\alpha]}_{n,\boldsymbol{\ell}}$ are given by $$\begin{aligned}
 P^{[\alpha]}_{0,\boldsymbol{0}}(\mathbf{x})=\,1,\label{eq:28}\quad P^{[\alpha]}_{1,\mathbf{e}_{j}}(\mathbf{x})=\,\frac{x_j}{\sqrt{\alpha}},\quad
 