First-order continuous- and discontinuous-Galerkin moment models for a linear kinetic equation: model derivation and realizability theory

TL;DR

This paper introduces new moment models for linear kinetic equations using finite element and discontinuous-Galerkin methods, demonstrating improved efficiency over classical models especially with non-smooth solutions.

Contribution

The paper develops two novel classes of moment models based on finite element and discontinuous-Galerkin methods, including realizability conditions and property analysis.

Findings

Models are more efficient than classical full-moment models in non-smooth cases.

Numerical tests confirm the effectiveness of the new models.

Realizability conditions are thoroughly investigated.

Abstract

We provide two new classes of moment models for linear kinetic equations in slab and three-dimensional geometry. They are based on classical finite elements and low-order discontinuous-Galerkin approximations on the unit sphere. We investigate their realizability conditions and other basic properties. Numerical tests show that these models are more efficient than classical full-moment models in a space-homogeneous test, when the analytical solution is not smooth.

Figures

Equations

$${\partial }_t\psi +\Omega \cdot {\nabla }_{\mathbf{x}}\psi +{\sigma }_a\psi ={\sigma }_s\mathcal{C}\left(\psi \right)+Q,$$

$$\mathcal{C}\left(\psi \right)=\underset{{\mathcal{S}}^2}{\int }K(\Omega ,{\Omega }^{\prime })\psi (t,\mathbf{x},{\Omega }^{\prime })\text{leavevmode}d{\Omega }^{\prime }-\underset{{\mathcal{S}}^2}{\int }K({\Omega }^{\prime },\Omega )\psi (t,\mathbf{x},\Omega )\text{leavevmode}d{\Omega }^{\prime }.$$

$$\psi (0,\mathbf{x},\Omega )$$

$$\psi (t,\mathbf{x},\Omega )$$

$$\Omega ={\left(\sqrt{1-{\mu }^2}\cos (\phi ),\sqrt{1-{\mu }^2}\sin (\phi ),\mu \right)}^T=:{\left({\Omega }_x,{\Omega }_y,{\Omega }_z\right)}^T,$$

$$\mathbf{u}=\underset{{\mathcal{S}}^2}{\int }\mathbf{b}\psi \text{leavevmode}d\Omega =:⟨\mathbf{b}\psi ⟩,$$

$$⟨\mathbf{b}{\partial }_t\psi ⟩+⟨\mathbf{b}{\nabla }_{\mathbf{x}}\cdot \Omega \psi ⟩+⟨\mathbf{b}{\sigma }_a\psi ⟩={\sigma }_s⟨\mathbf{b}\mathcal{C}\left(\psi \right)⟩+⟨\mathbf{b}Q⟩.$$

$${\partial }_t\mathbf{u}+⟨\mathbf{b}{\nabla }_{\mathbf{x}}\cdot \Omega \psi ⟩+{\sigma }_a\mathbf{u}={\sigma }_s⟨\mathbf{b}\mathcal{C}\left(\psi \right)⟩+⟨\mathbf{b}Q⟩.$$

$$\mathcal{H}(\psi )=⟨\eta (\psi )⟩$$

$$⟨\mathbf{b}\psi ⟩=\mathbf{u}.$$

$${\hat{\psi }}_{\mathbf{u}}=\underset{\psi :\eta (\psi )\in L_1}{\mathrm{argmin}}\left\{⟨\eta (\psi )⟩:⟨\mathbf{b}\psi ⟩=\mathbf{u}\right\}.$$

$$\mathbf{\alpha }(\mathbf{u}):=\underset{\tilde{\mathbf{\alpha }}\in {\mathbb{R}}^n}{\mathrm{argmin}}⟨{\eta }_{\ast }({\mathbf{b}}^T\tilde{\mathbf{\alpha }})⟩-{\mathbf{u}}^T\tilde{\mathbf{\alpha }},$$

$${\hat{\psi }}_{\mathbf{u}}={\eta }_{\ast }^{\prime }\left({\mathbf{b}}^T\mathbf{\alpha }(\mathbf{u})\right),$$

$${\partial }_t\mathbf{u}+{\partial }_x⟨{\Omega }_x\mathbf{b}{\hat{\psi }}_{\mathbf{u}}⟩+{\partial }_y⟨{\Omega }_y\mathbf{b}{\hat{\psi }}_{\mathbf{u}}⟩+{\partial }_z⟨{\Omega }_z\mathbf{b}{\hat{\psi }}_{\mathbf{u}}⟩+{\sigma }_a\mathbf{u}={\sigma }_s⟨\mathbf{b}\mathcal{C}\left({\hat{\psi }}_{\mathbf{u}}\right)⟩+⟨\mathbf{b}Q⟩.$$

$$\eta (\psi )=\psi \log (\psi )-\psi ,$$

$$\eta (\psi )=\psi \log (\psi )+\left(1-\psi \right)\log \left(1-\psi \right)$$

$$\eta (\psi )=\psi \log (\psi )-\left(1+\psi \right)\log \left(1+\psi \right)$$

$$\eta (\psi )=\frac{1}{2}{\psi }^2,$$

$${\hat{\psi }}_{\mathbf{u}}={\mathbf{b}}^T\mathbf{\alpha }(\mathbf{u})$$

$$\mathbf{\alpha }(\mathbf{u})={\mathbf{M}}^{-1}\mathbf{u}$$

$${\partial }_t\mathbf{u}+{\partial }_x{\mathbf{F}}_1\left(\mathbf{u}\right)+{\partial }_y{\mathbf{F}}_2\left(\mathbf{u}\right)+{\partial }_z{\mathbf{F}}_3\left(\mathbf{u}\right)=\mathbf{s}\left(\mathbf{u}\right),$$

$${\mathbf{F}}_1\left(\mathbf{u}\right)$$

$$\mathbf{s}\left(\mathbf{u}\right)$$

$${\partial }_t\psi +\mu {\partial }_z\psi +{\sigma }_a\psi ={\sigma }_s\mathcal{C}\left(\psi \right)+Q,t\in T,z\in X,\mu \in [-1,1].$$

$${\partial }_t\mathbf{u}+{\partial }_z⟨\mu \mathbf{b}{\hat{\psi }}_{\mathbf{u}}⟩+{\sigma }_a\mathbf{u}={\sigma }_s⟨\mathbf{b}\mathcal{C}\left({\hat{\psi }}_{\mathbf{u}}\right)⟩+⟨\mathbf{b}Q⟩.$$

$$\displaystyle\operatorname{span}({\mathbf{b}})=\big{\{}f\in L_{1}(V,\mathbb{R})\mid f\equiv\sum_{i=0}^{n-1}a_{i}b_{i},a_{i}\in\mathbb{R}\big{\}},$$

$${\mathbf{f}}_N={\left(1,\mu ,\dots ,{\mu }^N\right)}^T;$$

$${\mathbf{f}}_N={\left(P_0^0,P_1^0,P_2^0\dots ,P_N^0\right)}^T$$

$$h_i(\mu )=\{\begin{array}{ll}\text{mathbbm}{1}_{I_0}\frac{\mu -{\mu }_1}{{\mu }_0-{\mu }_1}& \text{ if }i=0,\\ \text{mathbbm}{1}_{I_{i-1}}\frac{\mu -{\mu }_{i-1}}{{\mu }_i-{\mu }_{i-1}}+\text{mathbbm}{1}_{I_i}\frac{\mu -{\mu }_{i+1}}{{\mu }_i-{\mu }_{i+1}}& \text{ if }0<i<k,\\ \text{mathbbm}{1}_{I_{k-1}}\frac{\mu -{\mu }_{k-1}}{{\mu }_k-{\mu }_{k-1}}& \text{ if }i=k,\end{array}$$

$${⟨{\hat{\psi }}_{\mathbf{u}}{h}_i⟩}_{I_{j-1}}$$

$${⟨{\hat{\psi }}_{\mathbf{u}}{h}_i⟩}_{I_j}$$

$$u_i$$

$${\mathbf{u}}_{I_j}={⟨{\mathbf{p}}_{I_j}\psi ⟩}_{I_j}={\int }_{I_j}{\mathbf{p}}_{I_j}\psi \text{leavevmode}d\mu .$$

Full text

††journal: arXiv.org††volume: Vi∩Vj

1 Introduction

Moment closures are a type of (non-linear) Galerkin approximation typically used in the context of kinetic transport equations. An infinite set of moment equations is defined by taking velocity- or phase-space averages of the kinetic density with respect to some basis of this space. A reduced description of the kinetic density is then achieved by truncating this hierarchy of equations at some finite order. Unfortunately, in most cases, the remaining equations are not closed and require information from the equations which were removed. The specification of this information, the so-called moment closure problem, distinguishes different types of moment models. In the context of linear radiative transport, the standard spectral method is commonly referred to as the $\text{P}_{N}$ closure [41], where $N$ is the degree of the highest-order moments in the model. It is basically a straight-forward Galerkin approximation on the unit sphere. The $\text{P}_{N}$ method is powerful and simple to implement, but does not take into account that the original function to be approximated, the kinetic density, must be non-negative. This often leads to physically meaningless solutions, as the $\text{P}_{N}$ solutions can, e.g., contain negative values for the local particle density.

In the context of radiative transport, entropy-based moment closures, the so-called $\text{M}_{N}$ models [44, 19], have all the properties one would desire in a moment method, namely positivity of the underlying kinetic density, hyperbolicity of the closed system of equations, and entropy dissipation [40]. Although there have been a lot of theoretical investigations about minimum-entropy models throughout the last fifty years, see, e.g., [38, 42, 12, 19, 30], practical implementations were limited to very low orders (e.g. $N=1$) [9, 46, 13, 14] as, generally, $\text{M}_{N}$ models are too expensive and more difficult to implement compared to, e.g., direct Monte Carlo or Discrete Ordinates simulations [25, 53, 5]. They require the (numerical) solution of a non-linear optimization problem at every point on the space-time grid, which tends to be very time-consuming. However, there has been renewed interest in these models recently due to their inherent parallelizability [27, 2].

The standard $\text{M}_{N}$ models use a polynomial basis on the velocity space. As a consequence, non-smooth kinetic densities (which very frequently occur in realistic problems) are poorly captured and often a very high moment order $N$ is needed for a sufficient approximation. To increase the accuracy of the $\text{M}_{N}$ models while maintaining the lower cost for small $N$, partial moment and mixed moment models have been developed [21, 20, 54, 47, 56]. They are based on a partition of the velocity space while keeping the moment order fixed, similar to some h-refinement for, e.g., finite element approximations [6].

While the partial moment models have been extensively studied for special cases (like half- or quarter-moments in one or two dimensions), we are unaware of any general investigation, especially in the fully three-dimensional setup.

In this paper, which is the first of two parts, we provide realizability theory for general first-order piece-wise discontinuous partial moments in one- and three-dimensional geometry, as well as their continuous analogue. An extensive numerical study using a second-order realizability-preserving scheme will be done in the second part. Here, the required realizability can cause problems for numerical methods, as standard high-order numerical solutions (in space and time) can destroy this property if not handled very carefully [62, 50, 4, 55, 52, 13, 46].

This paper is organized as follows. First, the transport equation and its moment approximations are given. Then, the new moment bases are presented and the available realizability theory is derived in one and three dimensions, followed by some numerical investigations of the approximation properties of our models. Finally, conclusions and an outlook on future work is given.

2 Modeling

We consider the linear transport equation

[TABLE]

The equation is supplemented with initial condition and Dirichlet boundary conditions:

[TABLE]

where $\mathbf{n}$ is the outward unit normal vector in $\mathbf{x}\in\partial X$.

Parameterizing $\Omega$ in spherical coordinates we obtain

[TABLE]

where $\varphi\in[0,2\pi]$ is the azimuthal and $\mu\in[-1,1]$ the cosine of the polar angle.

Definition 2.1**.**

The vector of functions ${\mathbf{b}}:\mathcal{S}^{2}\to\mathbb{R}^{n}$ consisting of $n$ basis functions $b_{i}$, $i=1,\ldots n$, of maximal order $N$ (in $\Omega$) is called an angular basis. The so-called moments $\mathbf{u}=\left(u_{0},\ldots,u_{n-1}\right)^{T}$ of a given density function $\psi$ are then defined by

[TABLE]

where the integration is performed component-wise.

Furthermore, the quantity $\rho=\rho(\mathbf{u}):=\left<\psi\right>$ is called the local particle density. Additionally, $\mathbf{u}_{\text{iso}}=\left<{\mathbf{b}}\right>$ is called the isotropic moment.

Equations for $\mathbf{u}$ can then be obtained by multiplying (2.1) with ${\mathbf{b}}$ and integration over $\mathcal{S}^{2}$, yielding

[TABLE]

Collecting known terms, and interchanging integration and differentiation where possible, the moment system has the form

[TABLE]

Depending on the choice of ${\mathbf{b}}$ the terms $\left<\Omega_{x}{\mathbf{b}}\psi\right>$, $\left<\Omega_{y}{\mathbf{b}}\psi\right>$, $\left<\Omega_{z}{\mathbf{b}}\psi\right>$, and in some cases even $\left<{\mathbf{b}}\mathcal{C}\left(\psi\right)\right>$, cannot be given explicitly in terms of $\mathbf{u}$. Therefore an ansatz $\hat{\psi}$ has to be made for $\psi$ closing the unknown terms. This is called the moment-closure problem.

In this paper the ansatz density $\hat{\psi}$ is reconstructed from the moments $\mathbf{u}$ by minimizing the entropy-functional

[TABLE]

The kinetic entropy density $\eta:\mathbb{R}\to\mathbb{R}$ is strictly convex and twice continuously differentiable and the minimum is simply taken over all functions $\psi=\psi(\Omega)$ such that $\mathcal{H}(\psi)$ is well defined. The obtained ansatz $\hat{\psi}=\hat{\psi}_{\mathbf{u}}$, solving this constrained optimization problem, is given by

[TABLE]

This problem, which must be solved over the space-time mesh, is typically solved through its strictly convex finite-dimensional dual,

[TABLE]

where ${\eta}_{*}$ is the Legendre dual of $\eta$. The first-order necessary conditions for the multipliers $\boldsymbol{\alpha}(\mathbf{u})$ show that the solution to (2.6) has the form

[TABLE]

where ${\eta}_{*}^{\prime}$ is the derivative of ${\eta}_{*}$.

Substituting $\psi$ in (2.4) with $\hat{\psi}_{\mathbf{u}}$ yields a closed system of equations for $\mathbf{u}$:

[TABLE]

This approach is called the minimum-entropy closure [39, 44, 43, 24]. The resulting model has many desirable properties: symmetric hyperbolicity, bounded eigenvalues of the directional flux Jacobian and the direct existence of an entropy-entropy flux pair (compare [39, 53]).

The kinetic entropy density $\eta$ can be chosen according to the physics being modeled.

As in [39, 27], we focus on Maxwell-Boltzmann entropy

[TABLE]

thus ${\eta}_{*}(p)={\eta}_{*}^{\prime}(p)=\exp(p)$. This entropy is used for non-interacting, classical particles as in an ideal gas. Other physically relevant entropies are [39]

[TABLE]

for particles satisfying Fermi-Dirac (e.g. fermions) (${\eta}_{*}(p)=\log(\exp(p)+1)$, ${\eta}_{*}^{\prime}(p)=\frac{\exp(p)}{1+\exp(p)}$) or

[TABLE]

for particles satisfying Bose-Einstein statistics (e.g. bosons) (${\eta}_{*}(p)=-\log(1-\exp(p))$, ${\eta}_{*}^{\prime}(p)=\frac{\exp(p)}{1-\exp(p)}$). Although photons are no classical particles but bosons, the approximation using the Maxwell-Boltzmann entropy is widely used due to its simplicity [44]. However, as in our special case the adaption to the Bose-Einstein case is rather straight-forward, we also include it in our demonstration code. Note that the modifications needed for the Fermi-Dirac case are more involved, as the additional requirement $\psi\leq 1$ changes many properties of the models. In particular, realizability (see below) has to be treated separately for this entropy.

Note that these physically relevant entropies all result in positive ansatz densities (2.8) (assuming ${\mathbf{b}}^{T}\boldsymbol{\alpha}(\mathbf{u})<0$ in the Bose-Einstein case). However, positivity is actually not gained for every entropy-based moment closure. Using the entropy

[TABLE]

the linear ansatz

[TABLE]

is obtained, leading to standard continuous/discontinuous-Galerkin approaches. In this case, the optimization problem can be solved analytically yielding

[TABLE]

where $M_{ij}=\left<b_{i}b_{j}\right>$ [53, 3]. If the angular basis is chosen as spherical harmonics of order $N$, (2.9) turns into the classical $\text{P}_{N}$ model [8, 10, 58], which thus can be considered a minimum-entropy model that does not guarantee positivity.

For convenience, we write (2.9) in the standard form of a non-linear system of hyperbolic balance laws:

[TABLE]

where

[TABLE]

To complete the definition of our moment method, we have to choose an angular basis. This will be done in the next section. Due to the notational complexity of the full three-dimensional setting we will first derive the models in slab geometry, which is a projection of the sphere onto the $z$-axis [58]. The transport equation under consideration then has the form

[TABLE]

The shorthand notation $\left<\cdot\right>=\int\limits_{-1}^{1}\cdot\leavevmode\nobreak\ d\mu$ then denotes integration over $[-1,1]$ instead of $\mathcal{S}^{2}{}$. Finally, the moment system is given by

[TABLE]

3 Angular bases in slab geometry

Definition 3.1**.**

Let ${\mathbf{b}}$ be an angular basis. The span of ${\mathbf{b}}$ is defined as

[TABLE]

where $V=[-1,1]$ in slab geometry and $V=\mathcal{S}^{2}$ in the full three-dimensional setting.

Definition 3.2**.**

Let ${\mathbf{b}_{1}}$ and ${\mathbf{b}_{2}}$ be two angular bases. We call them equivalent if $\operatorname{span}({\mathbf{b}_{1}})=\operatorname{span}({\mathbf{b}_{2}})$.

3.1 Full moments

The standard Galerkin approach, which results in the classical $\text{M}_{N}$ and $\text{P}_{N}$ models (using the entropies (2.10) and (2.11), respectively), chooses a polynomial basis ${\mathbf{b}}=\mathbf{f}_{N}$ on $[-1,1]$ as angular basis. There are two obvious options for $\mathbf{f}_{N}$, resulting in equivalent models:

• •

the monomial basis

[TABLE]

• •

the Legendre basis

[TABLE]

with the Legendre polynomials $P_{l}^{0}$, ${l}=0,\ldots,N$.

The monomial basis is attractive due to its simplicity, especially for realizability considerations [15]. On the other hand, the Legendre polynomials form an orthogonal basis of $L_{2}([-1,1],\mathbb{R})$ and are additionally eigenfunctions of the isotropic-scattering operator $\mathcal{C}$ if $K\equiv\frac{1}{2}$, diagonalizing the source-term. We thus use Legendre polynomials for the one-dimensional $\text{M}_{N}$ and $\text{P}_{N}$ models in our numerical experiments.

The full-moment approach is able to accurately approximate smooth density functions, often with just a few moments. However, it struggles with discontinuous densities which are abundant in practice. In addition, an obvious problem of full-moment models is that averaging over the complete angular domain may remove necessary information. In the following, we will thus investigate basis functions with local support in the angular domain. Although higher-order models are possible, we will stick to first-order approximation, i.e. piecewise linear basis functions.

3.2 Continuous piecewise-linear angular basis (Hat functions)

As mentioned above, the method of moments is nothing else than a Galerkin approximation of the kinetic equation (2.16) in the angular variable. The full-moment $\text{P}_{N}$ model is equivalent to the finite element method (FEM) with fixed angular mesh-size and increasing order (this method is often referred to as the $p$-version, related to the maximal order $p$ of the polynomial basis [6]). Another natural approach is the standard finite element method, where the polynomial order is fixed but the angular mesh is refined (often referred to as the $h$-version, where $h$ denotes the mesh-size [6]).

Let $\mathcal{I}_{k}=\mathcal{I}_{k}(\mu_{0},\ldots,\mu_{k})$ be the partition of $[-1,1]$ given by a set of $k+1$ angular “grid” points $-1=\mu_{0}<\mu_{1}<\ldots<\mu_{k-1}<\mu_{k}=1$ and corresponding intervals $I_{j}=[\mu_{j},\mu_{j+1}]$, $j=0,\ldots,k-1$. Given this partition, the continuous piecewise-linear basis functions $\mathbf{h}_{\mathcal{I}_{k}}=\left(h_{0},\ldots,h_{k}\right)^{T}$ (hat functions or B-splines of first order [18]) are defined as

[TABLE]

Definition 3.3**.**

The linear (using entropy (2.11)) and nonlinear model (using entropy (2.10)) with angular basis ${\mathbf{b}}=\mathbf{h}_{\mathcal{I}_{k}}$ will be called $\text{HFP}_{\mathcal{I}_{k}}$ and $\text{HFM}_{\mathcal{I}_{k}}$, respectively. If $\mathcal{I}_{k}=\mathcal{I}_{k}(-1,-1+\frac{2}{k},\ldots,-1+(k-1)\cdot\frac{2}{k},1)$ is the equidistant partition with $k$ intervals, we will refer to the basis as $\mathbf{h}_{n}$ and to the models as $\text{HFP}_{n}$ and $\text{HFM}_{n}$, respectively, where $n=k+1$ is the number of basis functions.

One of the most important properties of the hat functions is the following.

Lemma 3.4**.**

The hat functions form a partition of unity, i.e. $\sum\limits_{i=0}^{n-1}h_{i}\equiv 1$.

Corollary 3.5**.**

1. (a)

The full moment basis $\mathbf{f}_{1}$ is equivalent to $\mathbf{h}_{2}$. 2. (b)

The mixed moment basis **[22, 54]** $\mathbf{m}_{1}=\left(1,\mathbbm{1}_{[0,1]}\leavevmode\nobreak\ \mu,\mathbbm{1}_{[-1,0]}\leavevmode\nobreak\ \mu\right)^{T}$ is equivalent to $\mathbf{h}_{3}$.

Proof.

The claim follows immediately, observing that

1. (a)

$\mathbf{f}_{1}=\left(1,\mu\right)^{T}=\left(h_{1}+h_{2},-h_{1}+h_{2}\right)^{T}$, 2. (b)

$\mathbf{m}_{1}=\left(1,\mathbbm{1}_{[0,1]}\leavevmode\nobreak\ \mu,\mathbbm{1}_{[-1,0]}\leavevmode\nobreak\ \mu\right)^{T}=\left(h_{1}+h_{2}+h_{3},h_{3},-h_{1}\right)^{T}$. ∎

An advantage of the piecewise linear basis over higher order polynomials is that integrals including the minimum entropy ansatz density (using Maxwell-Boltzmann entropy) $\hat{\psi}_{\mathbf{u}}=\exp\left(\sum_{i}h_{i}\alpha_{i}\right)$ can be calculated analytically. Exemplarily, the $\text{HFM}_{n}$ integrals are given by

[TABLE]

Note that for $\alpha_{i-1}\approx\alpha_{i}$ or $\alpha_{i+1}\approx\alpha_{i}$ (which occurs for example in the isotropic state) the closed formulas become numerically instable. This can be avoided by replacing them in this situation with a Taylor expansion of suitable order.

3.3 Partial moments (discontinuous-Galerkin ansatz)

The idea of partial moments is not to average over the full domain at once but to partition the sphere (or its projection) into disjoint parts and define moments separately for every element of the partition [49, 21]. One model of this class, which has been successfully applied to radiative transfer in one dimension, is the half-moment approximation [20].

Let the partition $\mathcal{I}_{k}$ be defined as in Section 3.2. For every interval $I_{j}$ we define moments by

[TABLE]

The basis functions $\mathbf{p}_{I_{j}}$ are monomials on their corresponding intervals and zero elsewhere. Consequently,

[TABLE]

Note that evaluation on the boundary of $I_{j}$ is not defined. The set of basis functions is then given by the vector $\mathbf{p}_{\mathcal{I}_{k}}=\left(\mathbf{p}_{I_{0}}^{T},\ldots,\mathbf{p}_{I_{k-1}}^{T}\right)^{T}$.

Definition 3.6**.**

The linear (using entropy (2.11)) and nonlinear model (using entropy (2.10)) with angular basis ${\mathbf{b}}=\mathbf{p}_{\mathcal{I}_{k}}$ will be called $\text{PMP}_{\mathcal{I}_{k}}$ and $\text{PMM}_{\mathcal{I}_{k}}$, respectively. If $\mathcal{I}_{k}$ is the equidistant partition, we will refer to the basis as $\mathbf{p}_{n}$ and to the models as $\text{PMP}_{n}$ and $\text{PMM}_{n}$, respectively, where $n=2k$ is the number of basis functions.

Remark 3.7**.**

The $\text{PMM}_{2}$ model is equivalent to the full-moment $\text{M}_{1}$ model.

The advantage of considering only first-order moments is that, similar to the continuous piecewise-linear ansatz, the moment integrals in (2.5b) can be calculated exactly. Assume that the ansatz satisfies $\left.\hat{\psi}_{\mathbf{u}}\right|_{I_{j}}=\exp\left(\alpha_{j,0}+\alpha_{j,1}\mu\right)$. Then

[TABLE]

Remark 3.8**.**

We want to note that our bases satisfy the following relation:

[TABLE]

4 Angular bases in three dimensions

Like in one dimension, the full moment models are defined using a polynomial basis $\mathbf{f}_{N}$ over the whole velocity space $\mathcal{S}^{2}$. As a standard choice, we use real spherical harmonics of order $N$, resulting in $n=(N+1)^{2}$ moments for the $\text{P}_{N}$ and $\text{M}_{N}$ models.

To define (continuous and discontinuous) locally supported bases, we have to specify a partition of the domain. Albeit both approaches are not limited to this, we only consider moments on spherical triangles. A generalization to arbitrary convex spherical polygons is straightforward.

Let $\mathcal{T}_{h}$ be a spherical triangulation of $\mathcal{S}^{2}$ and \mathchoice{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}\in\mathcal{T}_{h} be a spherical triangle with vertices A_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}, B_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}} and C_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}} (or $A$, $B$, $C$ as short notation). Furthermore, let $K$ be the flat triangle spanned by the vertices $A$, $B$ and $C$, i.e.  \mathchoice{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}=g\left(K\right) with $g\left(\mathbf{x}\right)=\frac{\mathbf{x}}{\left\|\mathbf{x}\right\|_{2}}\leavevmode\nobreak\ .$ This is shown exemplarily in Figure 1.

In the following, we will use a dyadic refinement [7] of the quadrants/octants. This is achieved by subdividing every spherical triangle into four new ones, adding vertices at the midpoints of the triangle edges. This is shown in Figure 1 for one quadrant (black) and two refinements (red and green). Starting from the octants, after $r$ refinements, the triangulation $\mathcal{T}_{h}(r)$ consists of $\left|\mathcal{T}_{h}(r)\right|=2\cdot 4^{r+1}$ spherical triangles with $n_{v}(r)=4^{r+1}+2$ vertices (see Figure 2 for the case $r=1$).

Remark 4.1**.**

Previous results indicate that the dyadic refinement is better initialized using the vertices of an icosahedron instead of an octahedron [7]. We decided to use this variant regardless to include the quarter-moment model into our sequence of refinements. However, this is only for investigative reasons. Any practical application should start from the icosahedron.

4.1 Barycentric-coordinate basis functions

We consider basis functions defined using spherical barycentric coordinates [11, 35, 48]. These basis functions are associated with vertices of the triangulation $\mathcal{T}_{h}$.

Definition 4.2**.**

*Let the spherical triangle

¿

$\textstyle K$

be given. The functions h_{A_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}}, h_{B_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}} and h_{C_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}} are defined by*

h_{A_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}}\left(\Omega\right)+h_{B_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}}\left(\Omega\right)+h_{C_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}}\left(\Omega\right)=1* for every point \Omega\in\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}} (partition of unity).* 2. 2.

*At the vertices v\in\{A_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}},B_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}},C_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}\} of *

¿

$\textstyle K$

:

[TABLE]

and likewise for h_{B_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}} and h_{C_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}}, respectively (Lagrange property). 3. 3.

\Omega\in\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}* is the Riemannian center of mass with weights h_{A_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}}\left(\Omega\right), h_{B_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}}\left(\Omega\right), h_{C_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}}\left(\Omega\right) and nodes A_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}, B_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}, C_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}, respectively.*

Lemma 4.3** ([48]).**

1. (a)

The functions h_{A_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}}, h_{B_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}} and h_{C_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}} from Definition 4.2 exist uniquely. 2. (b)

For every interior point \Omega\in\operatorname{int}\left(\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}\right) and every v\in\{A_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}},B_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}},C_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}\} it holds that $h_{v}\left(\Omega\right)>0$ and $h_{v}\left(\Omega\right)\geq 0$ for every $\Omega\in\mathcal{S}^{2}$.

Details on how to evaluate the spherical barycentric coordinates are given in [48]. For the theoretical investigation of the resulting moment models we only need the properties stated above.

Numbering all vertices of $\mathcal{T}_{h}$ as $v_{0},\ldots,v_{n_{v}-1}$, the full set of basis functions is given as

[TABLE]

with h_{i}|_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}=h_{A_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}} if $v_{i}$ corresponds to vertex $A$ in

¿

$\textstyle K$

(and likewise for $B$, $C$) and h_{i}|_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}\equiv 0 on spherical triangles that do not contain $v_{i}$. Since the basis functions form a partition of unity, it follows that $\rho=\sum\limits_{i=0}^{n_{v}-1}u_{i}$. We show one example of such a basis function in Figure 2.

As above (see Definition 3.3), the corresponding moment models will be called $\text{HFP}_{\mathcal{T}_{h}}$ and $\text{HFM}_{\mathcal{T}_{h}}$, respectively. If $\mathcal{T}_{h}=\mathcal{T}_{h}(r)$ is the triangulation obtained by $r$ dyadic refinements of the octants, the models will be called $\text{HFP}_{n_{v}}$ and $\text{HFM}_{n_{v}}$, where $n_{v}=4^{r+1}+2$ is the number of vertices of the triangulation.

4.2 Partial moments on the unit sphere

Partial-moment models on the sphere have been introduced in reference [21] (although the authors restricted their investigation to quarter moments). We will only investigate first-order partial moments. Let

¿

$\textstyle K$

be a spherical triangle. We then define

[TABLE]

where \operatorname{int}\left(\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}\right) is the interior of

¿

$\textstyle K$

. As in one dimension, the value on the boundary of

¿

$\textstyle K$

remains unspecified. We will denote the moments with respect to \mathbf{p}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}} as \mathbf{u}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}=\left<\mathbf{p}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}\psi\right>. If \mathbf{u}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}={(u_{0},u_{1},u_{2},u_{3})}^{T} we define \mathbf{u}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|0\right|}=u_{0} and \mathbf{u}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|1\right|}={(u_{1},u_{2},u_{3})}^{T}.

The angular basis $\mathbf{p}_{\mathcal{T}_{h}}$ now consists of all element bases \mathbf{p}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}, \mathchoice{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}\in\mathcal{T}_{h}. The local particle density is then given by \rho=\sum\limits_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}\in\mathcal{T}_{h}}\mathbf{u}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|0\right|}. As above (see Definition 3.6), the corresponding moment models will be called $\text{PMP}_{\mathcal{T}_{h}}$ and $\text{PMM}_{\mathcal{T}_{h}}$, respectively. For $\mathcal{T}_{h}$= $\mathcal{T}_{h}$($r$), we will use the notation $\text{PMP}_{n}$ and $\text{PMM}_{n}$, where $n=4\cdot\left|\mathcal{T}_{h}(r)\right|=2\cdot 4^{r+2}$ is the number of moments.

Remark 4.4**.**

As in one dimension (see Remark 3.8), our bases satisfy the following relation:

[TABLE]

5 Realizability

Since the underlying kinetic density to be approximated is non-negative, a moment vector only makes sense physically if it can be associated with a non-negative density function. In this case, the moment vector is called realizable.

For higher-order numerical methods, realizability becomes a major challenge. Generally, the system of moment equations only evolves on the set of realizable moments. This property, though desirable since it is consistent with the original kinetic density, can cause problems for numerical methods. Standard high-order numerical solutions to the Euler equations, which indeed are an entropy-based moment closure, have been observed to have negative local densities and pressures [62]. This can be overcome by using a realizability-preserving limiter (see, e.g., [62, 4, 13, 46, 52]). However, this limiter requires exact knowledge of the set of realizable moments which often poses severe problems. In particular, for high moment orders in the full-moment context, calculating the limited approximation is numerically intractable [55].

In this section, we will investigate the structure of the realizable set and derive realizability conditions for our models, showing that they are particularly simple in all cases, which will later allow for an efficient high-order implementation in space and time for our minimum-entropy models $\text{HFM}_{n}$ and $\text{PMM}_{n}$.

5.1 Structure of the realizable sets

As already mentioned above, a moment vector is called realizable if it can be associated with a non-negative integrable density function. The set of realizable moments is called the realizable set.

Definition 5.1**.**

The (${\mathbf{b}}$-)realizable set $\mathcal{R}_{{\mathbf{b}}}$ is

[TABLE]

If $\mathbf{u}\in\mathcal{R}_{{\mathbf{b}}}$, then $\mathbf{u}$ is called ($\mathbf{b}$-)realizable. Any $\psi$ such that $\mathbf{u}=\left<{\mathbf{b}}\psi\right>$ is called a (${\mathbf{b}}$-)representing density.

Remark 5.2**.**

• (a)

The realizable set is a convex cone, and

• (b)

representing densities are not necessarily unique.

Unfortunately, the definition of the realizable set is not constructive, making it hard to check if a given moment vector is realizable or not. In the minimum-entropy context, realizability is further complicated by the fact that the optimization problem (2.6) only has a solution if the moment vector can be represented by a density of the form (2.8), i.e. if it lies in the ansatz set

[TABLE]

As the ansatz density ${\eta}_{*}^{\prime}\left({\mathbf{b}}^{T}\boldsymbol{\alpha}\right)$ is always strictly positive (assuming ${\mathbf{b}}^{T}\boldsymbol{\alpha}(\mathbf{u})<0$ in the Bose-Einstein case), we immediately get the inclusions

[TABLE]

In general, these inclusions might be strict and the ansatz set $\mathcal{A}_{{\mathbf{b}}}$ might be considerably smaller than the realizable set $\mathcal{R}_{{\mathbf{b}}}$. In that case, the minimum entropy ansatz would not be a reasonable approximation to the original kinetic density. Luckily, in the case of a bounded angular domain, it was shown in [30] (for Maxwell-Boltzmann entropy) that for some bases, including the full moment basis $\mathbf{f}_{N}$, the minimum entropy ansatz does not further restrict the set of realizable moments.

Theorem 5.3** ([30]).**

For the full moment basis we have

[TABLE]

The piecewise linear bases do not satisfy all of the properties required from the basis in [30]. More precisely, the assumption that

[TABLE]

is not valid. However, in the following we will show that this is not a problem as the ansatz set is still at least the interior of the realizable set.

Theorem 5.4**.**

Let ${\mathbf{b}}=\mathbf{h}_{n}$ or ${\mathbf{b}}=\mathbf{p}_{n}$. Then we have

[TABLE]

In particular, the ansatz set is dense in the realizable set, so each realizable moment vector can be accurately approximated by vectors in the ansatz set. However, note that numerical difficulties may occur when approaching the boundary of the realizable set [3].

Before proving Theorem 5.4, we will introduce a slightly more general realizable set based on distributions [57] instead of integrable functions which is often easier to work with for theoretical considerations. In the following, $V$ will always denote the angular domain ($[-1,1]$ in slab geometry and $\mathcal{S}^{2}$ in three dimension).

Definition 5.5**.**

A distribution on $V$ is a continuous linear functional mapping from the space $\mathcal{C}_{c}^{\infty}(V)$ of compactly supported smooth functions on $V$ to $\mathbb{R}$.

For a more precise definition and details on the topologies involved see, e.g., [29, 23]. In our case, $V$ is compact, so all functions on $V$ are compactly supported. Thus, for a smooth basis ${\mathbf{b}}$ (e.g. the full moment basis), evaluation of the distribution $\psi$ at the components of ${\mathbf{b}}$ is well-defined. The hatfunction and partial moment bases, however, are only piecewise smooth. We thus need the notion of a piecewise distribution.

Definition 5.6**.**

Let $\mathcal{P}=\{V_{0},\ldots,V_{k-1}\}$ be a partition of the angular domain $V$ with closed subsets $V_{j}$ such that $\bigcup\limits_{j=0}^{k-1}V_{j}=V$ and $=0\ \forall i,j$. A piecewise distribution on $\mathcal{P}$ is a collection $\psi=(\psi_{0},\ldots,\psi_{k-1})$ such that $\psi_{j}$ is a distribution on $V_{j}$, $j=0,\ldots,k-1$. We define the restriction of $\psi$ to $V_{j}$ as $\psi|_{V_{j}}\coloneqq\psi_{j}$.

In the following we will use $\mathcal{P}=\mathcal{I}_{k}$ and $\mathcal{P}=\mathcal{T}_{h}$ for the piecewise linear bases in one and three dimension, respectively, and $\mathcal{P}=\{V\}$ for the full moment basis.

We will abuse notation and write, for a distribution $\psi$ on $V$ and a smooth basis ${\mathbf{b}}$, the dual pairing as

[TABLE]

Note that if $\psi$ is an integrable function rather than a distribution, a distribution is defined by

[TABLE]

which shows that the notation $\left<\psi{\mathbf{b}}\right>$ for distributions is just a generalization of the notation for integrable functions. If $\psi$ is a piecewise distribution and $b$ is piecewise (with respect to $\mathcal{P}$) a smooth function on the angular domain, we further generalize the notation and define $\left<\psi b\right>\coloneqq\sum_{j=0}^{k-1}\psi|_{V_{j}}(b|_{V_{j}})$.

Remark 5.7**.**

Remember that we did not define the values of the partial moment bases on the boundary between two elements of the partition. If $\psi$ is an integrable function, these values can be chosen arbitrarily without changing the integrals as the boundary is a null set. When working with (piecewise) distributions, however, pointwise values have to be taken into account. We thus define the restriction $\mathbf{p}|_{V_{j}}$ of $\mathbf{p}$ to $V_{j}$ as the continuous extension of $\mathbf{p}|_{\operatorname{int}\left(V_{j}\right)}$ to $V_{j}$. In particular, $\mathbf{p}|_{V_{i}}$ and $\mathbf{p}|_{V_{j}}$ might differ on $V_{i}\cap V_{j}$.

Definition 5.8**.**

A distribution $\psi$ is called non-negative, written $\psi\geq 0$, if

[TABLE]

A piecewise distribution is called non-negative if all its restrictions to the elements of the associated partition are non-negative.

We are now ready to define the generalized realizable set.

Definition 5.9**.**

The generalized realizable set $\widetilde{\mathcal{R}}_{{\mathbf{b}}}$ is the set of moment vectors representable by non-negative (piecewise) distributions, i.e.,

[TABLE]

Note that by (5.3) it is clear that $\mathcal{R}_{{\mathbf{b}}}\subset\widetilde{\mathcal{R}}_{{\mathbf{b}}}$.

Example 5.10**.**

Consider the partial moment basis $\mathbf{p}_{4}={(\mathbbm{1}_{[-1,0]},\mu\mathbbm{1}_{[-1,0]},\mathbbm{1}_{[0,1]},\mu\mathbbm{1}_{[0,1]})}^{T}$ and the moment vector $\mathbf{u}={(2,0,1,0)}^{T}$. Then $\mathbf{u}=\left<\psi\mathbf{p}_{4}\right>$ with $\psi|_{[-1,0]}=2\delta$ and $\psi|_{[0,1]}=\delta$ and thus $\mathbf{u}\in\widetilde{\mathcal{R}}_{\mathbf{p}_{4}}$. However, $\mathbf{u}\not\in\mathcal{R}_{\mathbf{p}_{4}}$ because a non-negative integrable function $\psi$ with $\left<\psi\right>_{[-1,0]}=2$ will always yield $\left<\mu\psi\right>_{[-1,0]}<0$.

Remark 5.11**.**

For notational convenience, we assume that the Dirac distribution $\delta$ has full mass also on the boundary of integration.

If a (piecewise) distribution $\psi$ is (piecewise) a linear combination of Dirac deltas [26, 33], it is called atomic [15]. For the monomial basis in one dimension, every moment $\mathbf{u}\in\widetilde{\mathcal{R}}_{{\mathbf{b}}}$ realizable by a non-negative distribution can be represented by an atomic density [15]. Below, we will see that this is also true for our piecewise linear bases. We can use this to prove Theorem 5.4.

Proof of Theorem 5.4.

The proofs for the identity $\mathcal{A}_{{\mathbf{b}}}=\mathcal{R}_{{\mathbf{b}}}$ in [30], in particular the proof of [30, Lemma 6.2] using the null set assumption (5.2), also work for the piecewise linear bases if $\mathcal{R}_{{\mathbf{b}}}$ is replaced by $\mathcal{R}^{+}_{{\mathbf{b}}}{}$, which gives $\mathcal{A}_{{\mathbf{b}}}=\mathcal{R}^{+}_{{\mathbf{b}}}{}$. Further it was shown in [30] that

[TABLE]

(the proofs also hold for the piecewise linear bases). In particular, $\mathcal{R}^{+}_{{\mathbf{b}}}=\mathcal{A}_{{\mathbf{b}}}=\mathbf{u}(\mathbb{R}^{n})$ is open, so it remains to show that $\mathcal{R}^{+}_{{\mathbf{b}}}$ is dense in $\widetilde{\mathcal{R}}_{{\mathbf{b}}}$ (then, by the inclusion $\mathcal{R}_{{\mathbf{b}}}\subset\widetilde{\mathcal{R}}_{{\mathbf{b}}}$, $\mathcal{R}^{+}_{{\mathbf{b}}}$ is also dense in $\mathcal{R}_{{\mathbf{b}}}$). Let $\mathbf{u}\in\widetilde{\mathcal{R}}_{{\mathbf{b}}}$. Then it has an atomic representing density $\psi(\Omega)=\sum_{i}a_{i}\delta(\Omega-\Omega_{i})$ (see [15], Theorem 5.17, Theorem 5.22 and the proofs of Theorem 5.25 and Theorem 5.28 for the $\mathbf{f}_{N}$, $\mathbf{h}_{\mathcal{I}_{k}}$, $\mathbf{p}_{\mathcal{I}_{k}}$, $\mathbf{h}_{\mathcal{T}_{h}}$ and $\mathbf{p}_{\mathcal{T}_{h}}$ basis, respectively). Approximating each Dirac delta by a Dirac sequence of positive integrable functions $g_{i,j}$ converging to $\delta(\Omega-\Omega_{i})$ (in the sense of distributions), we obtain the sequence of positive densities $\psi_{j}=\sum_{i}a_{i}g_{i,j}$ whose corresponding moments $\mathbf{u}_{j}=\left<{\mathbf{b}}\psi_{j}\right>$ converge to $\mathbf{u}$. ∎

Remark 5.12**.**

The strictly positive Dirac sequences $g_{i,j}$ required in the proof of Theorem 5.4 can be obtained by choosing any of the compactly supported sequences given in [16], adding a constant $1/j$ on the domain of integration and renormalizing such that $\left<g_{i,j}\right>=1$ still holds. By [16, Theorem 1] (originally from [31]), the resulting sequence is again a Dirac sequence.

To summarize, we have the following structure of the realizable sets:

[TABLE]

where $\overline{\mathcal{R}_{{\mathbf{b}}}}$ is the closure of $\mathcal{R}_{{\mathbf{b}}}$. For the bases regarded here, the last inclusion also is an equality, as can be seen by the explicit descriptions provided in the following (see Theorems 5.13, 5.17, 5.22, 5.25 and 5.28).

Note that, in general, $\mathcal{R}_{{\mathbf{b}}}$ is neither open nor closed which makes it hard to give an explicit description of the realizable set. Instead, realizability conditions for its interior $\mathcal{R}^{+}_{{\mathbf{b}}}$ and its closure $\widetilde{\mathcal{R}}_{{\mathbf{b}}}$ will be given in the next subsections.

5.2 Slab geometry

Although completely solved, realizability conditions for the classical full-moment models turn out to be rather complicated [15, 36].

Theorem 5.13** ([15]).**

Let $\mathbf{u}\in\mathbb{R}^{N+1}$. Define the Hankel matrices

[TABLE]

and the rank of $\mathbf{u}$ with respect to the full-moment basis $\mathbf{f}_{N}$ as

[TABLE]

Using the definition $\mathbf{A}\geq\mathbf{B}$ if and only if $\mathbf{A}-\mathbf{B}$ is positive semi-definite, the following are equivalent.

$\mathbf{u}\in\widetilde{\mathcal{R}}_{\mathbf{f}_{N}}$** 2. 2.

$\begin{cases}\mathbf{A}(k)\geq\mathbf{0},\mathbf{A}(k-1)\geq\mathbf{C}(k)&\text{if }N=2k,\\ \mathbf{A}(k)\geq\mathbf{B}(k),\leavevmode\nobreak\ \mathbf{A}(k)\geq-\mathbf{B}(k)&\text{if }N=2k+1.\end{cases}$** 3. 3.

There exists an atomic representing density with $\mathrm{rank}_{\mathbf{f}_{N}}\left(\mathbf{u}\right)$ atoms.

“$1.\Leftrightarrow 2.$” also holds for the realizable set $\mathcal{R}_{\mathbf{f}_{N}}=\mathcal{R}^{+}_{\mathbf{f}_{N}}$ if “positive semi-definite” is replaced by “positive definite”.

Thus, to check whether a given moment vector is realizable with respect to $\mathbf{f}_{N}$, the definiteness of several Hankel matrices has to be tested. This is quite expensive numerically. In addition, a severe drawback of full-, partial- and mixed-moment, minimum-entropy models of more than first order is that the resulting integrals in the moment equations cannot be expressed in terms of elementary functions (not only in slab but also in the full geometry). This means that numerical quadrature is strictly necessary to solve these equations. Unfortunately, this has a strong impact on the realizable set and therefore also on the solution of (2.17) [2, 4], which further complicates testing realizability.

Fortunately, this does not hold for the first-order moment models where the integrals can always be evaluated. Moreover, even if the integrals are computed numerically, the realizable set does not change if a suitable quadrature is chosen (except for the partial moments in three dimensions, see Section 5.4). This allows to always use the analytical realizability conditions which will be derived in the following and which will turn out to be very easy to check (see Corollaries 5.19, 5.23, 5.27, 5.30).

5.2.1 Continuous piecewise-linear angular basis (Hat functions)

In order to derive a theorem similar to Theorem 5.13 for the hat function basis, we first have to introduce the rank of a moment vector with respect to this basis. We will start by introducing a decomposition of non-negative vectors into vector which have a single block of positive entries.

Definition 5.14**.**

Given a non-negative vector $\mathbf{u}\in{(\mathbb{R}_{\geq 0})}^{n}$, define the two index sets (start and end indices, respectively)

[TABLE]

where $u_{-1}=u_{n}\coloneqq 0$. Let $S(\mathbf{u})=\{s_{0},\ldots,s_{L-1}\}$ and $E(\mathbf{u})=\{e_{0},\ldots,e_{L-1}\}$ such that the indices $s_{l}$ and $e_{l}$ are sorted in ascending order. The positive blocks of $\mathbf{u}$ are then defined as

[TABLE]

for $l=0,\ldots,L-1$. We define the order of the positive block $\mathbf{u}_{l}^{+}$ as $\operatorname{ord}(\mathbf{u}_{l}^{+})=e_{l}-s_{l}+1$. The positive blocks give a unique decomposition

[TABLE]

of $\mathbf{u}$ (where the empty sum is defined as $\mathbf{0}$).

Example 5.15**.**

For $\mathbf{u}=(0,1,1,0,1,0,1,1)$ we have $S(\mathbf{u})=\{1,4,6\}$, $E(\mathbf{u})=\{2,4,7\}$, $L=3$, $\mathbf{u}_{0}^{+}=(0,1,1,0,0,0,0,0)$, $\mathbf{u}_{1}^{+}=(0,0,0,0,1,0,0,0)$ and $\mathbf{u}_{2}^{+}=(0,0,0,0,0,0,1,1)$ with $\operatorname{ord}{\mathbf{u}_{0}^{+}}=\operatorname{ord}{\mathbf{u}_{2}^{+}}=2$, $\operatorname{ord}{\mathbf{u}_{1}^{+}}=1$.

Based on this decomposition, we can now define the $\mathbf{h}_{\mathcal{I}_{k}}$-rank of a moment vector.

Definition 5.16**.**

Let $\mathbf{u}=\left(u_{0},\ldots,u_{k}\right)^{T}$ be given. If $u_{i}\geq 0$ for all $i=0,\ldots,k$ we define the rank of $\mathbf{u}$ with respect to the hat function basis by

[TABLE]

We are now ready to state the hat function basis analogue of Theorem 5.13.

Theorem 5.17**.**

Let $\mathbf{u}=\left(u_{0},\ldots,u_{k}\right)^{T}$ be given. Then the following are equivalent:

$\mathbf{u}\in\widetilde{\mathcal{R}}_{\mathbf{h}_{\mathcal{I}_{k}}}$** 2. 2.

$\mathbf{u}$* satisfies*

[TABLE] 3. 3.

There exists an atomic $\mathbf{h}_{\mathcal{I}_{k}}$-representing density with $r\coloneqq\mathrm{rank}_{\mathbf{h}_{\mathcal{I}_{k}}}\left(\mathbf{u}\right)$ atoms.

Proof.

• $1.\Rightarrow 2.$

Let $\mathbf{u}\in\widetilde{\mathcal{R}}_{\mathbf{h}_{\mathcal{I}_{k}}}$. Thus there exists a representing distribution $\psi\geq 0$. Since by construction $h_{i}\geq 0$, it holds that $\left<h_{i}\psi\right>\geq 0$. Therefore, (5.7) is necessary.

• $2.\Rightarrow 3.$

Let $n=k+1$. We first handle the case that $u_{i}>0$ for all $i$ and $n=2m$ is even. Thus, $\mathbf{u}=\mathbf{u}_{0}^{+}$ with $\operatorname{ord}(\mathbf{u}_{0}^{+})=n$. Define the atomic density

[TABLE]

It is easy to check that $\frac{u_{2i+1}\mu_{2i+1}+u_{2i}\mu_{2i}}{u_{2i+1}+u_{2i}}\in\left(\mu_{2i},\mu_{2i+1}\right)$. Additionally, $\left<h_{i}\psi\right>=u_{i}$, i.e. $\psi$ is a realizing distribution with $m$ atoms. The odd case $n=2m-1$ is a degenerated version of the even case (just by adding an interior node $\mu_{i}$, $0<i<k$, twice). This results in a non-uniqueness of the representing distribution (as the location of the additional node is free within the set of inner nodes).

Let now $\mathbf{u}=\sum_{l}\mathbf{u}_{l}^{+}$. If $\operatorname{ord}(\mathbf{u}_{l}^{+})=1$, a representing distribution with a single atom is given by $u_{s_{l}}\delta(\mu-\mu_{s_{l}})$. For $\operatorname{ord}(\mathbf{u}_{l}^{+})\geq 2$, we can get a representing distribution of the form (5.8) by ignoring the zero entries and treating $(u_{s_{l}},\ldots,u_{e_{l}})$ as a vector with strictly positive entries on the subpartition $\mathcal{I}_{\operatorname{ord}(\mathbf{u}_{l}^{+})-1}(\mu_{s_{l}},\ldots,\mu_{e_{l}})\subset\mathcal{I}_{k}$. A representing distribution for $\mathbf{u}$ is then given by adding the representing distributions of its positive blocks.

• $3.\Rightarrow 1.$

Trivial. ∎

Remark 5.18**.**

The definition of $\widetilde{\mathcal{R}}_{{\mathbf{b}}}$ uses piecewise distributions (see Definition 5.6), while we used distributions on the whole angular domain in the proof of Theorem 5.17. However, for atoms in the interior of an interval it is obvious how to get a piecewise distribution from the atomic distribution. Atoms on the interval boundaries can be arbitrarily assigned to one of the two distributions in the adjacent intervals due to the continuity of the hatfunction basis. In fact, since non-negative distributions are always of order 0 (see, e.g, [29, Theorem 2.1.7]) and thus can be evaluated on continuous functions, the definition of $\widetilde{\mathcal{R}}_{{\mathbf{b}}}$ by piecewise distributions is actually only required for the (discontinuous) partial moment bases.

Corollary 5.19**.**

$\mathbf{u}\in\mathcal{R}^{+}_{\mathbf{h}_{\mathcal{I}_{k}}}$* if and only if $u_{i}>0$ for all $i=0,\ldots,k$.*

Proof.

Follows directly from Theorem 5.17 using the the fact that $\mathcal{R}^{+}_{\mathbf{h}_{\mathcal{I}_{k}}}$ is the interior of $\widetilde{\mathcal{R}}_{\mathbf{h}_{\mathcal{I}_{k}}}$ (see 5.6). ∎

Remark 5.20**.**

Another $\mathbf{h}_{\mathcal{I}_{k}}$-representing distribution for $\mathbf{u}$ is

[TABLE]

5.2.2 First-order partial moments

For the partial moment basis, the rank is obtained by counting the intervals in the partition $\mathcal{I}_{k}$ where the representing distribution has non-zero mass.

Definition 5.21**.**

Let $\mathbf{u}_{\mathbf{p}_{\mathcal{I}_{k}}}=\left(u_{0,0},u_{1,0},\ldots,u_{0,k-1},u_{1,k-1}\right)^{T}$ be given. We define its rank with respect to the partial moment basis by

[TABLE]

The partial moment analogue of Theorem 5.13 is then easily obtained.

Theorem 5.22**.**

Let $\mathbf{u}=\left(u_{0,0},u_{1,0},\ldots,u_{0,k-1},u_{1,k-1}\right)^{T}\in\mathbb{R}^{2k}$ be given. Then the following are equivalent:

$\mathbf{u}\in\widetilde{\mathcal{R}}_{\mathbf{p}_{\mathcal{I}_{k}}}$. 2. 2.

$\mathbf{u}$* satisfies*

[TABLE] 3. 3.

There exists an atomic $\mathbf{p}_{\mathcal{I}_{k}}$-representing density with $r:=\mathrm{rank}_{\mathbf{p}_{\mathcal{I}_{k}}}\left(\mathbf{u}\right)$ atoms.

Proof.

As the partial moments are just full moments of order 1 on each interval, this follows immediately from the results in [15, Thm. 4.1+4.3]. ∎

Corollary 5.23**.**

$\mathbf{u}\in\mathcal{R}^{+}_{\mathbf{p}_{\mathcal{I}_{k}}}$* if and only if*

[TABLE]

Remark 5.24**.**

Possible $\mathbf{p}_{\mathcal{I}_{k}}$-representing densities for $\mathbf{u}\in\mathcal{R}^{+}_{\mathbf{p}_{\mathcal{I}_{k}}}$ include

[TABLE]

and

[TABLE]

5.3 Three dimensions

We will now turn to the fully three-dimensional setting. Since the rank of a moment vector with respect to the piecewise linear bases in three dimensions is much harder to define than in one dimension (among others, it depends non-trivially on the partition $\mathcal{T}_{h}$ of the sphere) we will restrict ourselves to deriving realizability conditions.

5.3.1 Spherical barycentric

Theorem 5.25**.**

Let $\mathbf{u}=\left(u_{0},\ldots,u_{n_{v}-1}\right)^{T}\in\mathbb{R}^{n_{v}}$ be a vector of moments. Then it is necessary and sufficient for the existence of a non-negative distribution $\psi$ which realizes $\mathbf{u}$ with respect to $\mathbf{h}_{\mathcal{T}_{h}}$ that

[TABLE]

Proof.

The proof works similar to the one-dimensional setting (compare Theorem 5.17). We first show that equation (5.12) is necessary. Let $\psi\geq 0$ be arbitrary. By Lemma 4.3(b) it holds that $h_{i}\geq 0$ which implies that $u_{i}=\left<h_{i}\psi\right>\geq 0$ and thus (5.12).

Assume that $\mathbf{u}$ satisfies (5.12). A representing distribution is given by

[TABLE]

Due to (4.1) it follows that

[TABLE]

∎

Remark 5.26**.**

We want to note that the above proof is not limited to barycentric coordinates but also works for other bases with values in $[0,1]$ and the Lagrange interpolation property (4.1).

Corollary 5.27**.**

$\mathbf{u}\in\mathcal{R}^{+}_{\mathbf{h}_{\mathcal{T}_{h}}}$* if and only if*

[TABLE]

5.3.2 Partial moments

Theorem 5.28**.**

Let \mathbf{u}=\left(\mathbf{u}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|0\right|},\mathbf{u}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|1\right|}\right)^{T}\in\mathbb{R}^{4} be a vector of moments. Then it is necessary and sufficient for the existence of a non-negative distribution $\psi\not\equiv 0$ which realizes $\mathbf{u}$ with respect to \mathbf{p}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}} that

[TABLE]

satisfies \boldsymbol{\phi}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|1\right|}\in\operatorname{conv}\left(\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}\right), where $\operatorname{conv}\left(\cdot\right)$ denotes the convex hull.

Proof.

The necessity of (5.14) follows immediately, since by assumption $\psi\geq 0$ is supported in

¿

$\textstyle K$

and integration is a linear operation.

To show the sufficiency of (5.14) we give a realizing distribution function under this assumption. The boundary of the convex hull \operatorname{conv}\left(\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}\right) is given by

¿

$\textstyle K$

and $K$. First, we give realizing distributions on these boundary parts. Due to the convexity of the convex hull and the linearity of the integral, every interior point can then be reproduced by a suitable convex combination of these two candidates.

Assume that \left\|\boldsymbol{\phi}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|1\right|}\right\|_{2}=1, i.e. \boldsymbol{\phi}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|1\right|}\in\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}. A realizing distribution is given by

[TABLE]

where $\boldsymbol{\delta}$ denotes the multi-dimensional Dirac-delta distribution. By assumption, \psi_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}} is supported in

¿

$\textstyle K$

(this is no longer true if \left\|\boldsymbol{\phi}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|1\right|}\right\|_{2}<1). Thus

[TABLE]

If \boldsymbol{\phi}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|1\right|}\in K (the opposite boundary of the convex hull), a realizing distribution is

[TABLE]

where $\lambda_{A}$, $\lambda_{B}$, $\lambda_{C}$ are the barycentric coordinates of \boldsymbol{\phi}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|1\right|} with respect to the vertices $A$, $B$ and $C$ on $K$, i.e. $\lambda_{A}+\lambda_{B}+\lambda_{C}=1$ and \lambda_{A}A+\lambda_{B}B+\lambda_{C}C=\boldsymbol{\phi}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|1\right|}. It immediately follows that $\psi_{K}$ is supported in

¿

$\textstyle K$

and realizes $\mathbf{u}$.

If \boldsymbol{\phi}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|1\right|} is in the interior of the convex hull, we find that \frac{\boldsymbol{\phi}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|1\right|}}{\left\|\boldsymbol{\phi}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|1\right|}\right\|_{2}} can be realized by (5.15) and the intersection point of the line from the origin $(0,0,0)$ and \frac{\boldsymbol{\phi}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|1\right|}}{\left\|\boldsymbol{\phi}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|1\right|}\right\|_{2}} with $K$ can be realized by (5.16). Certainly, \boldsymbol{\phi}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|1\right|} lies on the same line and can thus be realized by a convex combination of these two distribution functions. ∎

Remark 5.29**.**

Note that, similarly to the one-dimensional case, the $\mathbf{p}_{\mathcal{T}_{h}}$-representing non-negative piecewise distribution for a moment vector $\mathbf{u}$ might have zero mass on some of the spherical triangles. In this case \mathbf{u}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|0\right|}=0 and \mathbf{u}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}}^{\left|1\right|}=\mathbf{0}.

Corollary 5.30**.**

$\mathbf{u}\in\mathcal{R}^{+}_{\mathbf{p}_{\mathcal{T}_{h}}}$* if and only if*

[TABLE]

5.4 Numerically realizable set

In general, we cannot solve the integrals in (5.1) analytically and have to approximate them by a numerical quadrature $\mathcal{Q}$. For the full moments, this can have a severe impact on the realizable set [2, 4]. In the following, we will show that for the piecewise linear bases, except for the partial moments in three dimensions, these problems do not occur if a suitable quadrature is chosen.

Let $\mathcal{P}=(V_{0},\ldots,V_{k-1})$ be a partition of the angular domain $V$ in closed subsets as in Definition 5.6 and let $\mathcal{Q}=(\mathcal{Q}_{0},\ldots,\mathcal{Q}_{k-1})$ such that $\mathcal{Q}_{j}=\{(w_{i}^{j},\Omega_{i}^{j})\}_{i=0}^{{n_{\mathcal{Q}}}_{j}-1}$ is a quadrature on $V_{j}$, i.e., for an integrable function $b$

[TABLE]

is the approximation of the corresponding integral $\left<\cdot\right>$ with the quadrature rule $\mathcal{Q}$.

We define the numerically realizable set

[TABLE]

Note that pointwise values of ${\mathbf{b}}$ have to be available for (5.17) to be well-defined. For the partial moment bases, we thus use the continuous extension from $\operatorname{int}\left(V_{j}\right)$ to $V_{j}$ on each element of the partition (see Remark 5.7).

In general, $\mathcal{R}_{{\mathbf{b}}}^{\mathcal{Q}}$ is a strict subset of $\widetilde{\mathcal{R}}_{{\mathbf{b}}}$, i.e., there are some moments that are realizable analytically but cannot be realized with $\mathcal{Q}$. For the hat functions and one-dimensional partial moments, however, the numerically realizable set and the generalized realizable set agree for suitable quadratures.

Theorem 5.31**.**

Let ${\mathbf{b}}$ be an angular basis in one dimension with $n$ piece-wise linear (possibly discontinuous) basis functions, i.e., there exists a partition $\mathcal{I}_{k}(\mu_{0},\ldots,\mu_{k})$ of $[-1,1]$ (see Section 3.2) such that the restriction on the interval $I_{j}$ satisfies $\left.b_{i}\right|_{I_{{j}}}\in P^{1}(I_{{j}})$ for all $j=0,\ldots,k-1$ and $i=0,\ldots,n-1$. Additionally, let $\mu_{j}$ and $\mu_{j+1}$ be part of the node set of the numerical quadrature $\mathcal{Q}_{j}=\{(w_{i}^{j},\hat{\mu}_{i}^{j})\}_{i=0}^{{n_{\mathcal{Q}}}_{j}-1}$ on the interval $I_{j}$, i.e., $\mu_{j},\mu_{j+1}\in\{\hat{\mu}_{i}^{j}\leavevmode\nobreak\ |\leavevmode\nobreak\ i=0,\ldots,{n_{\mathcal{Q}}}_{j}-1\}$. Then we have

[TABLE]

Proof.

It is straight-forward to show that (generally) $\mathcal{R}_{{\mathbf{b}}}^{\mathcal{Q}}\subset\widetilde{\mathcal{R}}_{{\mathbf{b}}}$. Thus let $\mathbf{u}\in\widetilde{\mathcal{R}}_{{\mathbf{b}}}$. Since the basis is piecewise linear, we can find a representing density $\hat{\psi}$ such that $\hat{\psi}|_{I_{{j}}}(\mu)=\alpha_{0,j}\delta(\mu-\mu_{j})+\alpha_{1,j}\delta(\mu-\mu_{j+1})$ with suitably chosen $\alpha_{0,j},\alpha_{1,j}\in\mathbb{R}_{\geq 0}$ (compare (5.9) and (5.11)). Setting

[TABLE]

shows that $\mathbf{u}\in\mathcal{R}_{{\mathbf{b}}}^{\mathcal{Q}}$. ∎

With a very similar proof, it can be shown that $\mathcal{R}_{\mathbf{h}_{\mathcal{T}_{h}}}^{\mathcal{Q}}=\widetilde{\mathcal{R}}_{\mathbf{h}_{\mathcal{T}_{h}}}$ also in three dimensions if the quadrature contains the vertices $v_{i}$ of the triangulation.

For the partial moment basis in three dimensions, there is a numerical equivalent to the analytical realizability conditions (5.14):

Theorem 5.32**.**

Let \mathbf{u}=\left(\mathbf{u}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}_{j}}^{\left|0\right|},\mathbf{u}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}_{j}}^{\left|1\right|}\right)^{T}\in\mathbb{R}^{4}\setminus\{\mathbf{0}\} be a vector of moments, and let $\mathcal{Q}_{j}$ be a quadrature on \mathchoice{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}_{j}. Then \mathbf{u}\in\mathcal{R}_{\mathbf{p}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}_{j}}}^{\mathcal{Q}} if and only if \mathbf{u}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}_{j}}^{\left|0\right|}>0 and

[TABLE]

Proof.

In [3, Proposition 1] it was shown that the set of numerically realizable moments with density $\rho=1$ is the convex hull of the basis evaluations,

[TABLE]

Applying this to \frac{\mathbf{u}}{\mathbf{u}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{2.30208pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}_{j}}^{\left|0\right|}} and using the definition of \mathbf{p}_{\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.51207pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{3.2229pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}_{j}} gives the result. ∎

Note that Theorem 5.32 shows that the analytical and numerical realizable set differ for the three-dimensional partial moments as the quadrature would have to contain all points on

¿

$\textstyle K$

to reproduce the analytical realizability conditions.

6 Results

To test the approximation properties of the regarded models, we investigate convergence against some prescribed density function $\psi$.

For each basis ${\mathbf{b}}$, the moment vector $\mathbf{u}$ is calculated as $\mathbf{u}=\left<{\mathbf{b}}\psi\right>$. Then the ansatz densities are computed by (2.8). We use the backtracking Newton scheme from [3] (without change of basis). As an initial guess, we use the multipliers for the isotropic moment. Integrals are calculated using very fine quadratures to avoid numerical errors. In one dimension, for the $\text{M}_{N}$ models, the domain $[-1,1]$ is partitioned in $200$ equally spaced intervals and a Gauss-Lobatto quadrature of order 197 is used on each interval. The $\text{HFM}_{n}$ and $\text{PMM}_{n}$ models already come with the partition $\mathcal{I}_{k}$. We further subdivide $\mathcal{I}_{k}$ such that at least $200$ subintervals are available, then we use the same Gauss-Lobatto quadrature on each of the subintervals. In three dimension, we use a Fekete quadrature [60] of order $18$ on each spherical triangle on the most refined spherical triangulation for all models (including the $\text{M}_{N}$ models).

We follow the FAIR guiding principles for scientific research [61] and publish the code that generates the following results in [37]. The code allows to choose between Maxwell-Boltzmann and Bose-Einstein entropy. As the results with both entropies are qualitatively the same, we are only presenting the Maxwell-Boltzmann results in the following.

6.1 Slab geometry

In one dimension, we prescribe a Gaussian function

[TABLE]

as a smooth test case and a Heaviside function

[TABLE]

as a discontinuous test. Note that the minimum-entropy ansatz (2.8) is always strictly positive. The minimum entropy ansatz thus cannot represent density functions that are exactly zero in parts of the domain. Moreover, for very small (positive) values the multipliers $\boldsymbol{\alpha}$ tend to infinity (in absolute values) which makes it numerically infeasible to solve the dual optimization problem (2.7). We thus prescribe an isotropic “vacuum” density $\psi_{vac}$ (with local particle density $\rho=10^{-8}$) as a minimum instead of zero.

As a third test case, we use

[TABLE]

which describes two crossing beams of particles with velocities close to $-1$ and $0.5$, respectively. This is a hard test case for moment models, as it approaches the sum of two Dirac distributions. In particular, full-moment models of low order typically fail (e.g. $\text{M}_{1}$ produces the isotropic solution) in this specific setup [21, 54].

Exemplary ansatz functions are shown in Figures 4, 5 and 6. The convergence results can be found in Figure 3. Further plots for all tested models and tabulated errors are provided in the supplementary materials.

For the Gauss function, as expected [28], the $\text{P}_{N}$ models show exponential convergence whereas the $\text{M}_{N}$ models are exact (up to numeric errors) starting from order $2$. This is due to the fact that the Gaussian can be written as (2.8) if $N\geq 2$. For the piecewise linear models, second-order convergence is expected [17, 18], which is confirmed by the results both in $L_{1}$ norm and $L_{\infty}$ norm. The $\text{HFM}_{n}$ and $\text{PMM}_{n}$ approximations show similar shapes, dependent on the location of the ansatz intervals. For models with an odd number of intervals the peak of the Gauss function lies in the middle of an ansatz interval and is approximated by a constant value in that interval (due to its symmetry). These models thus show a plateau around zero. In contrast, models using an even number of intervals show a maximum at zero. In this test case, the linear $\text{HFP}_{n}$ and $\text{PMP}_{n}$ models do not differ much from their non-linear counterparts.

In the Heaviside test, the $\text{PMP}_{n}$ and $\text{PMM}_{n}$ models are exact if an even number of intervals is used (such that the jump occurs between two ansatz intervals). All other tested models and the $\text{PMP}_{n}$ and $\text{PMM}_{n}$ models for an odd number of intervals show the expected first-order convergence in $L_{1}$ norm [34, 59]. Notably, the $\text{HFM}_{n}$ models with an even number of intervals give errors that are several orders of magnitude lower than the other models though the rate of convergence is not improved. In $L_{\infty}$ norm, no convergence can be observed in this testcase, due to the well-known Gibbs phenomenon (see Figure 5). The $\text{M}_{N}$ and $\text{P}_{N}$ models both show strong oscillations around the reference solution. For higher orders, the frequency of the oscillations increases while their amplitude decreases, except at the discontinuity where the overshoot approximately stays the same. While the $\text{P}_{N}$ models are symmetric, the $\text{M}_{N}$ models mainly oscillate around the positive part of the Heaviside function due to their inherent positivity. The $\text{PMM}_{n}$ and $\text{PMP}_{n}$ models with an odd number of intervals are exact on all intervals except the one containing the discontinuity. In that interval, the $\text{PMM}_{n}$ models show a strong overshoot and the $\text{PMP}_{n}$ models show a symmetric under-/overshoot at both sides of the discontinuity. The $\text{HFM}_{n}$ and $\text{HFP}_{n}$ models show similar behavior, except that the under-/overshoots give rise to small oscillations in the other intervals due to the models’ continuity.

For the two crossing beams, approximation quality of the $\text{M}_{N}$ models initially rapidly increases with the moment order. $\text{M}_{1}$ and $\text{M}_{2}$ fail, giving (nearly) isotropic ansatz densities, while the $\text{M}_{3}$ solution starts to show features of the reference solution. $\text{M}_{4}$ is already quite close to the reference solution. For even higher orders, the approximation error decreases with second order in $n$. The first-order models struggle in this test case. Depending on the location of the ansatz intervals, the beam centered at $\mu=0.5$ is approximated either flattened out, sharpened or skewed to the left or right. The beam around $\mu=-1$ always is at the boundary of an interval and thus at least qualitatively met. Only for high moment numbers (starting at around $n=50$), the approximation quality significantly increases (with second order). Again, the linear $\text{HFP}_{n}$ and $\text{PMP}_{n}$ exhibit similar behavior. The $\text{P}_{N}$ models show significant oscillations and perform even worse than the piecewise first-order models for a large range of $n$. For very high orders, given that the crossing beams density regarded here is anisotropic but still smooth, the exponential convergence finally kicks in and the approximation errors rapidly decrease.

6.2 Three dimensions

In three dimensions, the tests again include a smooth Gaussian

[TABLE]

a discontinuous indicator function

[TABLE]

and a density function modelling two crossing beams of particles

[TABLE]

where $a=100$.

The convergence results are depicted in Figure 7. For completeness, we also show a selection of ansatz functions for Gauss (Figure 8), Square (Figure 9) and Crossing Beams test (Figure 10). Ansatz functions for all tested models and tabulated convergence results can be found in the supplementary materials.

The smooth Gauss function is reproduced exactly by all $\text{M}_{N}$ models. This is different to the slab geometry case since $\Omega_{x}^{2}+\Omega_{y}^{2}+\Omega_{z}^{2}=1$ and therefore this Gaussian on the unit sphere is actually just a first-order ansatz. Likewise, the $\text{PMM}_{n}$ models are exact. Due to quadrature errors and errors in solving the optimization problem, the actual computed error is slightly higher than numerical accuracy, which is especially visible for $\text{M}_{7}$ and $\text{M}_{8}$ ($64$ and $81$ moments, respectively). The $\text{HFM}_{n}$ ansatz is not able to model the Gaussian function exactly and results in a rather square-shaped ansatz density for the lower order models. The $\text{P}_{N}$ models show exponential convergence. For the piecewise linear models, quadratic convergence in the grid width is expected [28]. This corresponds to linear convergence in $n$, as the grid width approximately halves with each dyadic refinement, but $n=2\cdot 4^{k+2}$ after $k$ refinements for the partial moment basis and $n=2+4^{k+1}$ for the hatfunction basis. Both models converge with the expected order.

For the discontinuous square function, as in one dimension, first-order convergence in the grid width (for the piecewise linear models) is expected in the $L_{1}$-norm [28]. This corresponds to a convergence rate of $\frac{1}{2}$ in $n$. The expected convergence rate for the piecewise-linear models is achieved on average, with slightly varying values (from $0.3$ to $0.8$) between different refinement iterations. As in one dimension, on triangles containing or touching the discontinuity, artificial maxima are produced. On triangles away from the discontinuity, the exact solution is reproduced, with slight deviations in the case of the hatfunction models due to the continuity requirements. As it takes a very fine triangulation to resolve the discontinuity even the high-order models still show visible differences compared to the reference solution. The full moment models ($\text{M}_{N}$ and $\text{P}_{N}$) converge (on average) slightly slower, with an observed convergence rate of approximately $0.38$. Here, the low-order models cannot reproduce the square shape and are radially symmetric. Higher-order models closely mimic the square but show a regular oscillation pattern within the square domain. As in one dimension, the number of oscillations increases with the ansatz order. In $L_{\infty}$ norm, none of the models converge due to the famous Gibbs phenomenon.

6.3 Computational complexity

Solving the minimum entropy moment model (2.9) requires the (numerical) solution of the non-linear minimum entropy optimization problem (2.7) at every point on the space-time grid. Though inherently parallelizable (the optimization problems at a given time point are independent of each other) the computational cost often is prohibitively high for practical applications. Moreover, the numerical scheme has to guarantee that the moment vectors stay in the realizable set at all times as the optimization problem is not solvable otherwise. Given the complicated realizability conditions for the $\text{M}_{N}$ models (see Theorem 5.13), this is no easy task. Finally, the integrals needed during the optimization usually cannot be evaluated analytically for the $\text{M}_{N}$ models. Hence, a numerical quadrature has to be used which further complicates realizability.

The piece-wise linear models avoid many of these problems. The local support allows to use sparse or blocked containers in the implementation of the optimization algorithm. The realizability conditions are much simpler and, apart from the partial moment models in three dimensions, the numerically realizable set agrees with the analytically realizable set (see Section 5) which means that using a quadrature does not pose additional problems. In addition, analytical formulas for the integrals are available.

A thorough numerical investigation of the piece-wise linear minimum entropy moment models (2.9) will be done in a follow-up to this paper. Here, we only analyse the performance of solving a single optimization problem. We focus on the two crossing beams test cases which are representative of all our numerical tests. In general, the optimization algorithm converges faster for all models in the Gauss test cases and slower in the discontinuous ones but the relative behavior stays the same.

The results can be found in Figure 11. For the $\text{M}_{N}$ models, the computation times increase rapidly with higher order. The increase is quadratic in the number of moments $n$ in three dimensions and approximately of order $1.5$ in one dimension. In contrast, the execution times are basically independent of $n$ for the $\text{HFM}_{n}$ and $\text{PMM}_{n}$ models which are thus several orders of magnitude faster for large moment numbers. Here, due to the local support, only a fixed number of basis functions is non-zero at each point and thus the computational effort to, e.g., calculate the hessian matrix of the objective function only depends on the number of quadrature points. The zig-zag pattern in Figure 11a) results from our choice of quadrature intervals (see the beginning of Section 6) which does not yield exactly the same number of quadrature points for each model. A practical implementation would probably choose the quadrature suitable to the model, thus using a coarser quadrature for models with fewer moments. In that case, an additional linear increase in computational time with the number of quadrature points would be expected for all models.

7 Conclusions and outlook

We have derived minimum-entropy models for continuous and discontinuous piece-wise linear basis functions in one dimension and on the unit sphere. Additionally, we gave the corresponding realizability conditions, which also provide the set of moments for which these minimum-entropy models are hyperbolic and well-defined. Finally, numerical tests show the expected second- and first-order convergence to smooth and discontinuous kinetic densities, respectively.

For smooth tests, the standard minimum-entropy $\text{M}_{N}$ models show much better approximation properties than the new $\text{HFM}_{n}$ and $\text{PMM}_{n}$ models at equal moment number $n$. They should thus be preferred if smooth kinetic densities are expected. However, in practical applications, it is rarely known in advance that the kinetic densities are smooth, and non-smooth densities frequently occur. In that case, the new models give competitive approximations and are several orders of magnitude faster to compute. If in addition the ansatz triangulation is fitted to the discontinuities, the piece-wise linear models are even able to outperform the full moment models also with respect to approximation quality. Thus, approaches to adaptively choose the ansatz triangulation will be investigated further.

Comparing the two new model classes, continuous and discontinuous piece-wise linear basis functions generally perform very similar, though a direct comparison in three dimensions unfortunately is impossible due to the difference in the degrees of freedom (e.g. $\text{PMM}_{32}$ often yields a better result than $\text{HFM}_{18}$ but is outperformed by the next refinement $\text{HFM}_{66}$).

The sequel of this paper will deal with higher-order realizability-preserving numerical schemes for this type of moment models and an extensive numerical analysis for several classical benchmark problems. Additionally, we investigate the effects of numerical quadrature on the realizability set.

In future work, generalizations of the continuous piece-wise linear basis function to higher-order splines on the unit sphere [1] could be considered, as well as higher-order partial moments. However, realizability theory even for second order is challenging and has yet to be found, except for very few special cases like the full moment basis [32]. Additionally, the development of Kershaw closures as in [32, 45, 51, 52] for these new bases would be particularly interesting, as they provide a realizable moment method while avoiding the costly non-linear optimization problem of the minimum-entropy model.