Representation Theory Based Algorithm to Compute Boltzmann’s Bilinear Collision Operator in the Irreducible Spectral Burnett Ansatz Efficiently

Numerically solving the Boltzmann equation is computationally expensive in part due to the number of variables the distribution function depends upon. Another contributor to the complexity of the Boltzmann Equation is the quadratic collision operator describing changes in the distribution function due to colliding particle pairs. Solving it as efficiently as possible has been a topic of recent research, e.g. Cai and Torrilhon (Phys Fluids 31(12):126105, 2019. https://doi.org/10.1063/1.5127114), Wang and Cai (J Comput Phys 397:108815, 2019. https://doi.org/10.1016/j.jcp.2019.07.014), Cai et al. (Comput Fluids 200:104456, 2020. https://doi.org/10.1016/j.compfluid.2020.104456). In this paper we exploit results from representation theory to find a very efficient algorithm both in terms of memory and computational time for the evaluation of the quadratic collision operator. With this novel approach we are also able to provide a meaningful interpretation of its structure.


Introduction
Solving the Boltzmann equation numerically is an area of ongoing research. There are already existing results based on Spectral-Fourier Approach by Pareschi and Russo [1] and Gamba and Tharkabhushanam [2], with recent work improving some of its shortcomings such as [3] on the preservation of conservations or [4] on the stability of such methods. Another ansatz is based on the spectral Hermite ansatz introduced by Grad (1949) [5] with early work [6] and more recent work such as [7,8] [9] able to calculate the bilinear collision operator based on the spectral Hermite ansatz. However, their method is computationally expensive. In response, Cai et al. [10] developed and implemented a more efficient algorithm based on the spectral Burnett ansatz. This ansatz has also been looked into analytically [11] and numerically [12]. Inspired by Struchtrup's (2005) tensorial approach to the collision operator in Equation 6.35 in [13], we chose a different approach for calculating the Spectral-Burnett approximation of the Boltzmann collision operator. The main difference lies in the exploitation of features of the basis set, which is adapted to the irreducible subspaces with respect to the orthogonal group of the polynomial space as it consists of real solid spherical harmonics multiplied with Laguerre polynomials. While we also implemented the mostly analytical calculation of collision coefficients for a variety of potentials, the focus of this paper is the algorithm and corresponding numerical code, which uses these coefficients to solve the space-homogeneous Boltzmann equation for any distribution. Utilizing the properties of the irreducible subspaces, compared to the previous works we are able to significantly reduce the memory and computation time required for calculating the numerical solution. We could achieve this by virtue of the uniqueness up to a constant of linear maps between the irreducible subspaces. Representation theory allowed us to isolate each bilinear map, the constant for each map that encodes particle collisions and identify all bilinear maps, that evaluate to zero due to the mathematical properties of the irreducible subspaces in question. Our method provides a deep understanding of the underlying structure of the bilinear collision operator, allows for a very flexible setup for numerical computations and therefore is ideal for a purposeful reduction of computational effort.
The structure of this paper is as follows: The next section provides an overview of the collision operator and the potentials used in this paper. Sections 3 and 4 describe the application of representation theory to the collision operator, while Sect. 5 explains the representation theoretically inspired algorithm and its implementation. In Sect. 6 we apply our algorithm and also compare the results with [9]. For an overview of the applied representation theory refer to the "Appendix".

Boltzmann's Collision Operator
To describe gas flows, the Boltzmann equation is known to be valid for a wide range of situations. In particular, it describes non-equilibrium gas flows accurately for Knudsen numbers in the transition regime between the range of Navier-Stokes-Fourier equations and the kinetic regime. The Boltzmann equation describes the evolution of the gas distribution function f , which in general depends on seven variables-time t, position space x and velocity space c.
As done in previous papers we look at the space-homogenous Boltzmann equation to be able to focus solely on the collision operator: where f gives the time-dependent velocity distribution of particles, and S( f , f ) is the bilinear collision operator, see e.g. the textbooks [14,15], modelling the binary collisions of the particles in the gas: where n is the unit vector orthogonal to the collisional plane, the relative velocity g = c − c 1 , the average velocity h = 1 2 (c + c 1 ), and the primed velocities c , c 1 , g = c − c 1 with ||g || = ||g||, h = h are the velocities after the collision. Additionally, we assume the initial average velocity of the gas to be 0, and as this is one of the conserved moments it will always be 0. The solution of (1) converges to a Gaussian centered around the average velocity. Our assumption will allow us to exactly represent this solution with our numerical ansatz.
B(|g|, χ) = |g| d 2 4 with d a constant for the diameter of the colliding particles. (5) • Its generalization given by the variational hard spheres (VHS) model, • The inverse power law potentials (IPL), i.e. with η > 3 determining the power law exponent (7) and W 1 the positive root of We use the change of variables given in [7] to evaluate W 0 (χ).
We identify two types of assumed potentials corresponding to whether or not B depends on the angle χ. Additionally, η can be chosen from (3, ∞). Note that the IPL potential models can be derived from assuming a specific potential between the colliding particles, whereas the VHS model only allows that for η = ∞, in which case the two models correspond to each other. In this sense, the power potentials can be seen as more realistic models and the VHS models as a simplified version.
In this paper we will work with different examples of these models and compare them. In particular, we work with the classical χ-dependent Maxwell Potential with η = 5, the χ-independent and hence simplified Maxwell potential, as well as the hard power potential with η = 10. As we will explore in Sect. 4.2.2, these choices reflect different levels of computational difficulty both in a priori calculations as well as in the numerical evaluation.

Expansion of the Distribution Function
To solve the space homogeneous Boltzmann equation numerically, we follow [16] to expand the distribution function with an orthonormal basis set ϒ α (c), where α is typically a multi-index in three indices determining the basis polynomial α , f (eq) (c) is the global equilibrium or Maxwellian function and the scalar product is chosen for the weighted function space W eq : L 2 (R 3 , R, 1 f (eq) dc), see e.g. [16]: to allow us to identify the weights w α as moments of the distribution function. Here, we shift the velocity-coordinates such that the medium velocity is zero and rescale the variables so that the temperature is 1, giving us Generalizing this approach, the distribution function f (t, c) is approximated as a timedependent polynomial in c of degree M ∈ N or less, multiplied by the global Maxwellian, Gaussian-like equilibrium distribution , implicitly choosing the approximation space for the distribution function f (t, c) ∈ V eq with where R[c x , c y , c z ] is the space of all real polynomials in three variables and N = 1 6 (M + 1)(M + 2)(M + 3) the dimension of P M . We can equip the polynomial space with a scalar product corresponding to the setup (10):

Expansion of the Collision Operator
Inserting the expansion of the distribution function (9) into the space homogeneous Boltzmann equation (1) and using the linearity of the bilinear collision operator yields: We now project this expression onto V eq spanned by orthonormal basis functions ϒ α . The left hand side leads to while the right hand side evaluates to Overall, we find that the rate of change of the distribution function's weights depends on the collision vector s α , The collision vector in turn is determined by the bilinear operation (15) between the weights w α of the distribution function and a basis-dependent collision tensor

Structure of the Approximation Space
The remaining choice are the N basis polynomials α (c), with which we span P M , which has great influence on the structure and sparsity of the collision tensor and hence on the required computational time and memory to obtain s α . For example Grad [5] or more recently Cai [9] choose Hermite polynomials in three coordinates, since these are known to be orthogonal with respect to the corresponding scalar product ·, · P . However, structure and sparsity of the resulting collision tensor can still be improved.
We want to use a basis that exploits a key feature of the collision operator: It describes a physical process in the three dimensional velocity space, requiring it to be, like all physical processes, compatible with isometric transformations of the coordinate system, i.e. with O 3 the group of orthogonal matrices from R 3×3 : While this identity is often accepted, we included a proof in the "Appendix" for completeness, as this statement is central for this paper. This restriction allows us to employ results from the representation theory of O 3 and reduce the bilinear maps with which we express s α to only those bilinear maps between the O 3 -irreducible subspaces of P M .

Decomposition of P M into Its Irreducible Subspaces
With the space of homogeneous and harmonic polynomials of degree n ∈ { 0, 1, · · · , M } and dimension 2n + 1 we have the building blocks for the irreducible subspaces of P M [17]. To make it clear when an index describes the degree of a harmonic polynomial, we call n the anisotropic index.
We also define the one-dimensional space of homogeneous isotropic polynomials of second degree and denote with c 2 its element where λ = 1. Ignoring the chosen scalar product from Eq. (13) for now allows us to decompose P M into a direct sum of homogeneous polynomials with respect to some corresponding scalar product, We can check that this decomposition is complete by counting dimensions: Figure 1 sorts these subspaces column-wise according to their homogeneous degree deg = 2a + n and row-wise according to their anisotropic index n. Subspaces with same anisotropic index n are isomorphic to each other and the decomposition is not unique within each set of isomorphic subspaces. Subspaces with different anisotropic indices are non-isomorphic to each other and therefore necessarily pairwise orthogonal with respect to any O 3 invariant scalar product, see Remark 9.5 in this paper's "Appendix". Now we come back to the chosen scalar product in (13). The isomorphic subspaces shown in Fig. 1 are not orthogonal with respect to ·, · P . We therefore orthogonalize each set of isomorphic subspaces with a linear combination of suitable powers of (c 2 ) determined by associated Laguerre polynomials: where a H n are the orthogonalized irreducible subspaces. We call a the isotropic index and are left with the basis choice for H n .

Basis Choice Adapted to the Irreducible Subspaces
Here we employ the real-valued solid spherical harmonics Y l n : R 3 → R of n-th degree with the directional index l ∈ { −n, −n + 1, · · · , n − 1, n }, related to the real spherical harmonics Y l n : Note, that we can draw parallels to indices in quantum mechanics where the spherical harmonics are complex, and the principal, azimuthal and magnetic quantum numbers have similar mathematical roles to the isotropic index a, anisotropic index n and directional index l, respectively. This leads to the following choice for α : with (a, n, l) ∈ N 0 × N 0 ×{ −n, · · · , n } specifying the multi-index α, and the normalization coefficient given by

Discussion of the Basis Set and Irreducibility
Instead of immediately choosing a basis set, we first looked at the irreducible subspaces of P M with respect to O 3 and collected them as spaces of homogeneous polynomials in Fig. 1. Upon finding that we have multiple isomorphic subspaces, all row-wise collected in Fig. 1, we linearly combined them via Laguerre polynomials to obtain a H n , the orthogonalized irreducible subspaces. Then we chose basis polynomials for each a H n , but again using the isomorphism between spaces within one row, it was enough to choose a set of basis polynomials for each lowest degree irreducible subspace per row, i.e. for the spaces H n of harmonic polynomials of degree n.
We have chosen real solid spherical harmonics, since they are easily accessible, orthogonal and adapted to the chain of subgroups O 3 ≥ O 2 ≥ O 1 with c z the special symmetry direction, i.e. for an O 2 -invariant problem rotated so that c z is the rotation axis for O 2 , all coefficients with l = 0 are zero.
The basis set containing all basis polynomials with 2a + n ≤ M spans the full space P M . Note, that only basis polynomials with n = 0 are rotationally invariant. However, similar to the decomposition of a natural number into its prime factors, this basis is adapted to the decomposition of P M into its smallest invariant subspaces. Here "smallest" is defined as an invariant subspace U , that cannot be further decomposed into invariant subspaces other than U itself and the zero space, similar to a prime number only being devisable by itself and 1. These subspaces are called irreducible. Figure 1 shows the decomposition of P M into its homogeneous irreducible subspaces, while we use orthogonalized subspaces with respect to ·, · P given for a pair of chosen isotropic index a and anisotropic index n by with a dimension of 2n + 1. Note that we have recovered the spectral Burnett ansatz with real spherical harmonics. First results with the spectral Burnett basis have been developed previously in [10][11][12]18], but we are the first to exploit the underlying structure of the collision operator in our implementation. We also implement the nearly fully analytic calculation of the tensor coefficients and carry out numerical tests for large M.

Decomposing the Collision Operator
Equation (15) is the orthogonal projection of S( f , f ) onto the approximation space V eq . Equivalently we can identify the bilinear map C : P M × P M → P M with coefficients given by the collision tensor in terms of the basis functions of P M , now defined as C αβγ = C anl,bnl,cñl := anl , S( bnl f (eq) , cñl f (eq) ) P , i.e. with p, q ∈ P M and their respective expansion coefficients w anl , u anl we can express the map C via C( p, q) = anl bnl cñl C anl,bnl,cñl w bnl u cñl anl .
Not only is C a bilinear map, it inherits the property of the collision operator S that it is linear with respect to O 3 , compare to Eq. (18): Such a map is called natural. This restriction has been studied in representation theory and implies a decomposition of C into a sum over isotropic and anisotropic indices C : P M × P M → P M = abc nnñ (29) with S nnñ abc ∈ R some constants that, given a set of maps nnñ abc , distinguish between different maps C. Each nnñ abc : P M × P M → P M is a non-zero natural bilinear map given by with orthogonal projections P (n,a) : P M → a H n , isomorphisms ψ (n,a) : a H n → H n , embedding E (n,a) : a H n → P M , and the universal natural bilinear map This map Q nnñ is either zero or unique up to a constant, depending on the combination of anisotropic indices: Q nnñ is zero unless n +ñ − n ∈ 2N 0 , and 0 ≤n +ñ − n ≤ 2 min(n,ñ).
Note, that a decomposition of a (multi-) linear map between vector spaces into a sum of maps between subspaces as done in Eq. (29) by itself is generally possible. Equation (28) is crucial for the decomposition of C in (29), requiring nnñ abc and Q nnñ to be natural as well, and ultimately restricting all possible bilinear and natural maps Hn × Hñ → H n to be the same up to a constant with the additional condition (32).

Decoding the Structure of the Collision Tensor
Overall, the previous result implies the decomposition of the collision tensor into two tensors: mirroring the structure of irreducible subspaces of P M . Q nnñ lll and S abc nnñ have two distinct roles for the overall collision tensor. The coupling tensor Q nnñ lll collects the coefficients of the map Q nnñ . It is independent of the potential and encodes the purely mathematical coupling of Yln and Ylñ to Y l n . The impact tensor S abc nnñ encodes the physics of the particle collisions.

The Coupling Tensor
The coupling tensor Q nnñ lll encodes the natural bilinear map, unique up to a constant factor, from two harmonic polynomial spaces of degreesn andñ to a third harmonic polynomial space of degree n, see Eq. (31). Crucially, two natural bilinear maps between different, but isomorphic irreducible subspaces are identical up to a constant factor and hence all can be expressed in terms of one Q nnñ . Referring to Fig. 1 and knowing that all subspaces of the same row are isomorphic to each other, we can also picture Q nnñ as the bilinear map between three rows. This universality is the key in our algorithm, allowing us to achieve high efficiency both in terms of computational time and memory usage. We can calculate the values of the coupling tensor for example with The coupling tensor is related to the Clebsch-Gordon-Coefficients [19]. While the Clebsch-Gordon-Coefficients are more common in quantum physics, where they are usually used for three harmonic spaces of complex polynomials, the coupling tensor is the equivalent for real harmonic polynomials. It shares the property that it is zero for most combinations of anisotropic and directional indices (n,n,ñ, l,l,l), but does not have as many restrictions as its complex counterpart. In particular, the conditions in Eq. (32) hold for Q nnñ lll , and if they are met, then Q nnñ 000 is unequal to zero. Keeping in mind the uniqueness up to a constant of the coupling tensor, we choose one instance of Q nnñ for each triplet (n,n,ñ), calculate the values of the coupling tensor once and hard-code them in the implementation as function .

The Impact Tensor
In contrast to Q nnñ lll , the impact tensor S abc nnñ contains all information about the physics of the collision. While from the mathematical point of view the constant factor distinguishing two maps Q nnñ with same anisotropic indices (n,n,ñ) may seem irrelevant, it is clear that in the physical application these factors are the main interest. Here the conservation rules come into play as zeroes for particular combinations of anisotropic and isotropic indices a, b, c, (n,n,ñ), but this is also where different potentials lead to different values in the Table 1 Overview of a few choices for η and χ -dependence in B(|g|, χ) and their corresponding potential models Realistic power potential VHS model Hard power potential Hard VHS model η = ∞ Hard spheres potential Note, that for η = ∞ the power potential corresponds to the VHS model, and therefore the Hard spheres potential is χ -independent coefficients of the impact tensor. Notice, that the impact tensor is independent of the directional indices (l,l,l), so with our chosen map Q nnñ we can compute its coefficients with for combinations of anisotropic indices (n,n,ñ) fulfilling the conditions in (32). Depending on the potential, we can solve the integration in (35) fully or almost fully analytically with a computer algebra system, see the git-repository [20] for further details.
Referring to Table 1, we can now see where each chosen potential falls on the spectrum of high and low computational difficulty of both a priori evaluation, i.e. the calculation of the impact tensor itself, and runtime evaluation of the collision vector.
Note, that there is some relation between these two types of difficulty. The a priori evaluation only determines a set of numbers that make up the impact tensor, whereas the runtime evaluation mostly depends on the number of non-zero impact tensor entries and not on the specific non-zero entry S abc nnñ . As already noted by Grad [5], for the potentials with η = 5 the additional condition 2a applies, which stems from the orthogonality of the basis functions and the independence of the kernel from c and c 1 . This significantly reduces the amount of non-zero impact tensor coefficients, leading to comparatively low computational difficulty both a priori and during runtime. The a priori computational difficulty of the impact tensor coefficients additionally is determined by the χ-dependency of the kernel B. All χ-independent models are of lower difficulty, with the simplified Maxwell model the lowest difficulty due to (36). On the other side of the spectrum lies the hard power potential with η = 10, since it is χ-dependent and condition (36) does not apply, requiring the evaluation of the full set of S abc nnñ . This potential also has high runtime computational difficulty, since almost all impact tensor coefficients are non-zero.

Full Generalization of the Collision Operator
Notice, that the decomposition of the collision tensor into the impact and coupling tensor not only implies a particular implementation to calculate the collision vector, but it also gives us a unique understanding of the role of each coefficient. In particular, we can now generalize the collision operator to any symmetric, bilinear map between real polynomials in three variables with physical meaning, i.e. invariance under transformation of the coordinate system: Such a map would not necessarily follow the conservation laws but still follow the decomposition into coupling and impact tensor. To obtain an upper limit for the computational cost of calculating the collision vector, we can assign the impact tensor random numbers for all combinations of anisotropic and isotropic indices that fulfil Eq. (32). When testing the performance with this setup, we shall call such a potential general.

Implementation
The task is to implement the tensor-vector multiplication between the expansion coefficients of the collision tensor in the form of S abc nnñ Q nnñ lll and the expansion coefficients w anl of the distribution function. The result is saved as a vector with expansion coefficients s anl , so we want to find an efficient algorithm to compute the collision vector s, i.e. (37)

Custom Implementation Adapted to the Given Tensor Structure
Looking at the calculation task in (37), we can separate it into operations that need to be carried out for all combinations of anisotropic and isotropic indices. Beginning by setting the collision vector s to zero, we can rewrite the task as (38) Note, how the impact tensor is now outside of the sum, only remaining as a constant factor.

Implementing the Coupling Tensor
In fact, we can now isolate the map between the irreducible subspaces of different anisotropic and isotropic indices in terms of an operation that only involves the coupling tensors and the directional indices: (39) Here the anisotropic indices are only relevant to indicate the degree of the harmonic subspaces and hence the length of the involved vectors. Note, that we changed the positions of the indices and font of the coefficients of the distribution function and collision vector to indicate that this is an operation on a small subsection of the overall vectors.
Apart from the constant factor determined by the impact tensor, this operation is independent of the isotropic indices (a, b, c), so we can give it as a general operation between three generic vectors ς n , un, vñ of lengths 2n + 1, 2n + 1, 2ñ + 1, respectively: We obtain Eq. (39) by replacing β with the corresponding value of S abc nnñ . Since Q nnñ lll are constants, the above operation (40) can be hard-coded without loss of generality of the overall implementation. We do so by giving each triple of anisotropic indices (n,n,ñ) a unique case numer q and saving the expansion coefficients of both the distribution function w anl and the collision vector s anl in a specific order, so that we seperate it into coherent 2n + 1 long sections, or sub-vectors, w an . We also order the w anl such that the directional index l starts at −n and ends at n, see Fig. 2. This reflects the splitting of the polynomial space into its irreducible subspaces: Each section w an contains the coordinates of w in one particular irreducible subspace characterized by the anisotropic and isotropic indices a and n. The position of each section w an is given by a corresponding pointer p an , which locates the first of the 2n + 1 entries. The collision vector is ordered similarly.
Furthermore, we know that the following symmetry relations must hold: Therefore, we used Mathematica to generate the C++ function that implements Eq. (40) for all combinations of anisotropic indices (n,n,ñ) withn ≤ñ. Input for this method is β as double , q as integer and the three vectors ξ n , ζn, ςñ as three double pointers , , . The method switches between all relevant triplets of anisotropic indices (n,n,ñ), where q indicates the triplet and is thus used in the case distinction. For each case C++-code corresponding to the operation in (40) was generated. Here the three pointers are used as starting positions for the three vectors.
Note, that this code does not check, whether the three vectors actually have been allocated, all input is assumed to be valid in order for the computation to be as time-efficient as possible. Furthermore, it is always assumed thatn ≤ñ, the input vectors should be ordered accordingly (Fig. 3).

Implementing the Impact Tensor
In contrast to Q, the coefficients S abc nnñ can vary depending on the potential. We therefore give the impact tensor in a text file as input for the code. Each line of this text file is an input list containing the necessary information for , compare to Eq. (40) and see Fig. 4: The average timings are given for one timestep, i.e. they contain 4 evaluations of the collision operator. They have been calculated by measuring the calculation time of 1000 timesteps with a generic input distribution. Notice, that we are able to provide an upper bound for the computational requirements with the general potential that has non-zero entries for every impact tensor coefficient. Note, that the coupling tensor is stored as compiled C++ code and requires for M = 20 roughly 38 MB To calculate the collision vector, i.e. to compute (38), the C++-code iterates through each line of the impact tensor file and applies with the case q, β = S abc nnñ and the three vectors determined by p an , p bn , p cñ . Utilizing (41), we can shorten the input list by stating twice the value of S abc nnñ when (b,n) = (c,ñ), leaving the line withn <ñ and omitting the corresponding line of S acb nñn .

Discussion of the Implementation
Effectively, we have now implemented the map nnñ abc from Eq. (30): Structuring the array of moments w anl into subsections corresponding to the irreducible subspaces and using the positions given by the pointers p bn and p cñ as starting positions for the operation corresponds to φ (n,b) • P (n,b) and φ (ñ,c) • P (ñ,c) .
itself is the implementation of the map Q nnñ , only that it operates on the section of the vector s anl indicated by the pointer p an , implementing E (n,a) • φ −1 (n,a) . This setup is very flexible. It is easily adapted to different potentials, theories and initial conditions, i.e. the initial occupancy of w anl of the distribution function. Omitting lines belonging to irrelevant combinations of (a, b, c, n,n,ñ) for the given setup, we can build calculation-specific input files for the collision operator. In Table 2 we have listed the required memory space for the impact tensor and the average calculation time of one time step for a selection of different potentials and theories. All computation is done sequentially on a CPU model Intel(R) Core(TM) i7-8650U. We note that the hard power potential with η = 10 is very close to the general potential in terms of computational requirements. In Table 3 we have listed the number of coupling coefficients for different M and the corresponding sizes of the compiled program. We note, that the memory requirement for the compiled program is much larger than that for the text file of the impact tensor.

Comparison to an Algorithm for General Tensor-Vector Multiplications and Current State of the Art
Given that the ultimate goal of this work is to produce an implementation of the evaluation of S( f , f ) within a numerical code, that is as fast as possible, the question is, whether the above described implementation is better than an algorithm developed specifically for the computation of tensor-vector multiplications. The efficient implementation of tensor The memory requirements of the coupling tensor are given by the corresponding sizes of the compiled C++ code operations is a field of on-going research in high performance computing and currently there are many different options available.
To check how the algorithm we developed compares to other algorithms, we chose the Tensor Algebra Compiler (TACO) [21] as a case study, as it provides code generation in C++ and seems to perform similarly well to other libraries. Here we found that TACO performed best when given the full collision tensor C αβγ with every triple of indices (a, n, l) mapped to one counter index. We found the computational times comparable but slightly slower than our implementation. We also compare memory requirements for the full tensor C αβγ with the irreducible Burnett approach and found that for a general potential, the full tensor C αβγ is much more demanding than what is needed for the impact tensor and together. However, when measuring the memory requirements for the Maxwell potential, we initially found that the memory requirements for the full tensor C αβγ are less than the requirements for . This is because saving as compiled program is more memory intensive than a file with only numbers, and because we still compared to the general coupling tensor which has many more cases than actually necessary-for M = 20 out of the 1331 possible cases only 505 are actually required for the Maxwell potential, leading to a significantly smaller program of only 6.6MB. This reduction is due to the additional condition on the indices (36).
Assessing the performance of the irreducible Burnett algorithm would be incomplete without looking at the performance of the Burnett spectral method in [10] obtained from a single thread computation on a CPU model Intel Xeon E5-2697A V4. In Table 4 we can compare the three algorithms for the Maxwell potential. Additionally, for the hard Maxwell potential with η = 10 and M = M 0 = 20, we find from Figure 15 in [10], that one evaluation of the collision operator takes about 90ms, whereas our algorithm takes about 50ms. Note however, that a different computational setup with different CPU models was used for each computation. We therefore conclude, that the computational time of both methods is similar for the IPL potential with η = 10, while the memory cost is greatly reduced with the irreducible Burnett method, which requires about 51MB, whereas the spectral Burnett method needs 784MB for just the sparse Maxwell potential. Table 4 also allows us to compare performance for the Maxwell potential of both methods and here we can see a significant difference both in computational time and memory usage.
Overall, we are satisfied to conclude that our implementation utilizes the underlying structure of the collision operator well enough to outperform the current state of the art. The time is given for M = 20, 100 iterations and initial function given by the BKW solution for the Spectral Burnett algorithm and a generic input distribution for the other algorithms. Note, that all algorithms use precomputed coefficients for the collision operator and run sequentially. In these timings the irreducible Burnett algorithm has not been simplified to the homogeneous case, see Table 2 for these specs. The timings and memory usage for the spectral Burnett algorithm are taken from [10] with computation done on a single thread with CPU model Intel Xeon E5-2697A V4, while the other computations ran sequentially on a CPU model Intel(R) Core(TM) i7-8650U

Comparison to the Spectral Burnett Method in [10]
The spectral Burnett method used in the ansatz in [10] is based on complex solid harmonics.
Here the presented algorithm exploits the relation for the directional indices-the collision tensor's entries can only be non-zero if l =l +l holds. By contrast, the irreducible Burnett ansatz uses real solid harmonics, for which this condition does not apply. Instead, we find a sparse coupling tensor with band structures, since the real solid harmonics are a simple linear combination of their complex counterparts. Crucially, the main difference between the two approaches lies in the decomposition of the collision tensor into the impact and coupling tensor. The algorithm of the irreducible Burnett ansatz exploits both this decomposition and the sparsity of the coupling tensor. It is possible, that the performance can still be improved by applying the decomposition to the basis with complex spherical harmonics.

Comparison to the Spectral Fourier Approach
The irreducible and the spectral Burnett ansatz share a similar approach to discretize the velocity space through moments. This method is based on picking a finite polynomial degree and using a finite-dimensional approximation space for the velocity-dependency of the distribution function. In contrast, the spectral Fourier approach discretizes the velocity space by choosing some finite 3-dimensional velocity domain and representing that with a grid of equally distant points. The collision operator is a weighted convolution. [2] and [10] compared its computational costs to the presented algorithm and found comparable performance. Fixing the velocity domain introduces problems in the form of aliasing effects, choice of cut off and conservation properties. On the other hand, moment methods have difficulty in approximating some distributions, an example for this will be demonstrated in Sect. 6.3. Notably, the spectral Fourier approach is connected to the representation theory of abelian groups, such as O 2 or the group of translations. In fact, the functions e 2πiθ k used for the expansion each span a one-dimensional irreducible subspace [22].

Using the Irreducible Burnett Ansatz
As we will see in the numerical application, not all initial conditions are computable in the Burnett basis. However, moment equations allow a relatively low number of numerical degrees of freedom to describe physically relevant situations. Moments as variables encode physical quantities, usually those of low order are the most interesting, allowing us to interpret a subset of the numbers representing the distribution function in a meaningful way. Moment equations come with a level of complexity, and the irreducible Burnett ansatz is no exception. Its implementation is challenging compared to other numerical methods for the Boltzmann equation such as direct simulation monte carlo methods [23] or the aforementioned spectral Fourier approaches. This is due to the low level C++ code of the coupling tensor, which needs to be generated using scripts or similar for a high polynomial degree M due to the large number of different cases. We generated the function with Mathematica, and while the generating code can be checked, the overall code of cannot be verified by a human, only tested for plausibility in a small number of cases. However, this paper must be seen in the context of moment approximations for the Boltzmann equation, see e.g. [13,24,25]. In this context, via one way or the other, the collision operator will always be implemented as some matrix-vector product in the linear case, or a tensor-vector-vector product in the bilinear case. While the implementation of such a (bi-) linear operation can be verified by a human, for large degree M the matrix or tensor entries themselves similarly cannot be checked due to their number.
Solving just the spatially homogeneous Boltzmann equation is incomplete, the overall goal is to provide a method, with which the collision operator can be evaluated efficiently. However, the application of the irreducible Burnett ansatz to the inhomogeneous Boltzmann equation is out of the scope of this paper and we refer to [26] for further details.

Numerical Application
The irreducible Burnett ansatz is the real pendant to the spectral Burnett ansatz, additionally mathematically decomposing the collision tensor. While the spectral Burnett ansatz in [10] uses complex solid harmonics, the distribution function is real for their method as well. Therefore we expect its properties to be inherited from the spectral Burnett ansatz, including the rate of convergence. All differences between results of both methods must stem from purely numerical effects resulting from computational accuracy. This section aims to ensure the correct implementation of the method. To do so, we will roughly estimate spectral convergence by comparing lower orders of approximation to the numerical solution with M = 20 for all cases and times. Additionally we compare to the analytical solution where possible.
Let f denote the exact distribution,f the numerical solution with degree M = 20 and momentsŵ,f any other numerical solution with smaller degree and momentsw, f := f −f and their difference. Using the L 2 -norm implied by the scalar product on V eq from Eq. (10), f is simply given by the standard vector norm of the difference between their moments, Unlike the exact solution f , which we can only access in the first test case, f is available for all times. It is also a good indicator for spectral convergence, as we find for the residual Note that for t = 0,f is the result of the orthogonal projection of f and therefore f −f is in the orthogonal complement of V eq , implying that f ⊥ ( f −f ). Therefore, at t = 0 holds, so when f decreases, so does f −f . We will calculate f −f at t = 0 to find the difference's order of magnitude, but show plots of f for both t = 0 and t = 0, as we are mostly interested in the evolution of the coefficients calculated by the irreducible Burnett ansatz.
We consider three initial distributions to demonstrate the range of numerical application of our method: • BKW-Solution, where the initial function can be visualized as a 3-dimensional hollow sphere. • Bi-Gaussian distribution, where the initial function consists of two overlapping Gaussians. • Discontinuous distribution, where the initial function consists of two different states that are seperated at a flat plane in velocity space and that haven't interacted with each other yet.
We will look both at the effect of the choice of potential, as well as the difference between various numerical approximations, mainly given by the choice of different maximal polynomial degrees M. The equilibrium solution to the Boltzmann equation, the Gaussian distribution, lies in the chosen approximation space V eq for any M. With the our basis the Gaussian distribution corresponds to the state w 000 = 1 and all other moments equal to zero. Note, that all presented initial distributions evolve towards the Gaussian distribution with time for all presented potentials. Therefore, for large times we expect no difference between any two solutions, regardless of their initial distribution. This also implies that the residual between any two solutions decreases with increasing time. As a result, the focus should be on early time scales, because differences show up in early times and decrease with large time. Additionally, conservation of mass, momentum and energy imply constraints on the first coefficients, making them irrelevant for further consideration:

BKW Solution
The BKW solution is one of the few analytical solutions to the homogeneous Boltzmann Equation with Maxwell-potential and has been originally developed by Krupp (1967); Bobylev et al. (1976) [27]. As done in [9], we use the dimensionless formulation from Ernst [27]: where  The initial distribution f (0, c) can be envisioned as a 3D hollow sphere. In particular, it is obvious, that this distribution is rotationally invariant, allowing us to plot it in cylindrical coordinates without loosing information. Since our basis is adapted to O 3 , we can reduce the degrees of freedom to only those coefficients that correspond to rotationally invariant basis functions, i.e. only weights w anl with (a, n, l) = (a, 0, 0) can be non-zero (Fig. 5). Combined with the restraints from Eq. (42), this allows us to take a look at the remaining nine time-dependent moments, all of higher order with w 200 the lowest one, to judge how well the analytical and numerical solutions match without loosing any information.
Given that Cai et al. observed very good correspondence in [9,10], and that the approximations differ in basis choice at most, we expect to obtain similar results. Looking at Fig. 6, we observe that this expectation is met. We can also see spectral convergence of our method . We observe faster than exponential convergence, which is demonstrated by the downward sloping curves on the semi-log scale. However, there is no reason to expect that this behaviour translates from the very special BKW test case to the approximation of a general distribution.
Note, that for various M the difference between the evolution of the lowest order timedependent moment w 200 , which all compared theories share, is in the order of 10 −14 . This is due to the restriction in Eq. (36) imposed by the Maxwell potential, which for this particular case simplifies to With a = 2, the only possible contributions can come from (b, c) ∈ { (0, 2), (1, 1) }, which are already taken into account by the lowest order ansatz with M = 4. A corresponding argument applies to any higher moment-in this specific test case the increase of polynomial degree cannot have an influence on the lower moments, but it does add accuracy to the overall distribution, as shown in Fig. 7.
While the difference in basis choice and understanding of the underlying representation theory seems benign, in this particular case it means we have reduced the roughly 1700 perceived degrees of freedom to just 9. This also affects the number of relevant impact coefficients: For a rotationally invariant distribution and M = 20 we only need 34 impact coefficients S abc nnñ and hence the evaluation time for one time step is only a few ms. Figure 8 compares the evolution of the moments for different potentials. By omitting any line in the impact tensor file where b = 0, we generated the linear collision operator. In Fig. 8a we can clearly see a difference between the linear Maxwell potential and the analytical solution, in particular in the higher moments. We also calculated the evolution of the coefficients for the simplified Maxwell model and found that it corresponds well to the  Figure 8b compares the residual between the moments obtained from the numerical solution with M = 20 and the moments given by projecting the analytical solution onto the approximation space V eq with M = 20 for the three different potentials. We can see that both the Maxwell and the simplified Maxwell potential produce very accurate moments with error in the order of 10 −13 , while the linear potential differs significantly from the analytical solution. The plot also shows f −f , i.e. the difference between the analytical and numerical solution obtained with the Maxwell potential for the ansatz with M = 20 in the full function space equipped with the same scalar product as V eq . The residual starts at 8.9e−4 for t = 0 and we can see it steadily decreasing with time, as the solutions converge to the same Gaussian.
We conclude, that for rotationally invariant cases, the simplified Maxwell model produces similar results as the exact Maxwell model, and linearising the collision operator introduces Fig. 9 Bi-Gaussian distribution: Exact and with M = 20 approximated inital distribution function together with two views of their difference in cylindric coordinates a noticeable error. We also see even in this simple example the limit of spectral Burnett ansatz-moments can be accurately calculated, but the actual distribution function can only be approximated with some error depending on the degree M.
We also calculated the evolution of the coefficients for the simplified Maxwellian model and found that it corresponds to the analytical solution as well. The linearised simplified Maxwell model differs to the analytical solution similarly to the linearised Maxwell model. We conclude, that for rotationally invariant cases, the simplified Maxwell model is as good as the exact Maxwell model, and linearising the collision operator introduces a small but noticeable error (Fig. 9).

Bi-Gaussian Distribution
Again we take the initial distribution from [9] f (0, c) = 1 Due to the distribution's mirror symmetry w.r.t the x, y-plane, only even anisotropic indices n can result in non-zero projection coefficients in our basis. This is reflected in the projection of a monomial. Using the computer algebra system Mathematica, the projection of a monomial leads to a product from which it is easy to find a corresponding condition-any uneven degree leads to 0: 1 (a, b, c), with some additional function K 1 (a, b, c). We further see, that the distribution is We note, that with all the coefficients of uneven anisotropic indices being zero in the beginning and due to the conditions (32), we again may drop many lines from the impact tensor file, albeit not as many as in the BKW case. Figure 11 shows the residual between the ansatz with M = 20 and lower order approximations with the maximal polynomial degrees M ∈ { 10, 12, 14, 16, 18 }.
Here we again see spectral convergence of the solution, indicating that this test case can also be well approximated by the spectral Burnett basis. Figure 10 indicates that the smallest and largest orders M = 10 and M = 20 both have similar looking distributions, with their largest difference located at the tips and overlap of the two Gaussians. We can further see from both Figs. 10 and 11 that the difference between approximations decreases with time as discussed in the beginning of this chapter. f −f at t = 0 evaluates to 1.6e−3, the same order of magnitude as the first case. Overall, this indicates that the initial distribution is well approximated.
Lastly, we compare the simplified and actual Maxwell potential in Fig. 12, where we look at the evolution of 19 out of the 66 coefficients as examples to demonstrate the range of possibilities. We note that for some moments the Maxwell and simplified Maxwell potential correspond to each other well, while for other moments those calculated with the Maxwell potential converge faster to 0. This trend is observed for most of the not displayed coefficients as well. However, we also note that some coefficients grow first, before converging to zero.
Here we see, that the coefficients of the Maxwell potential remain closer to zero than those of the simplified Maxwell potential. Overall, we conclude that the simplified Maxwell potential can lead to a significantly different solution than the Maxwell potential in setups that are not rotationally invariant.

Discontinuous Distribution
Using the dimensionless formulation from [9] for this case we have Projection of a monomial, evaluated by Mathematica, again leads to a product from which we can infer that certain index combinations lead to zero: a, b, c), again with some additional function K 2 (a, b, c) However, finding exclusion rules for the projection coefficients is less obvious here. They can only be non-zero if n is uneven or zero. We again find that the distribution function is O 2 -symmetric, requiring the projection coefficients to be zero unless l = 0. Note in contrast to the bi-Gaussian case, that even though at the beginning the relevant coefficients have an uneven or zero anisotropic index, two uneven anisotropic indices can couple to an even anisotropic index, see (32). So the total set of relevant coefficients w anl for M = 20 is given by: isotropic index a ∈ {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10} ,   i.e. one coefficient per irreducible subspace, which makes 121 coefficients in total for M = 20. However, because all combinations of isotropic and anisotropic indices are possible, this case requires the full impact tensor and basis set, which has 1771 moments. Due to the O 2 -symmetry, most of them are zero and will always be zero for this test case. Looking at Fig. 13, we can clearly see that the discontinuity at t = 0 is not well approximated, even with a large degree of M = 20. This is reflected in f −f , which at t = 0 evaluates to 1.7e−1, a much larger value than in the other two cases.
Comparing the evolution of the distribution function in Fig. 14 for approximations of degree M = 10 and M = 20, we again find that the two approximations get closer to each other with time, as both converge to the same Gaussian. Figure 15 shows the residual between the ansatz with M = 20 and lower order approximations, each denoted by their maximal polynomial degree M ∈ { 10, 11,12,13,14,15,16,17,18 } .
Here we see a rather slow spectral convergence and that there is little gain of accuracy when increasing M from an uneven to an even degree, which is also why M = 19 was omitted. We propose that this is due to the initial condition, in which most of the coefficients with even 2a + n are zero. Together with the badly approximated discontinuity shown in Fig. 13, this shows one of the limits of the approximation space. This is the first test case where all lines of the impact tensor file are contributing to the overall evolution of the distribution function and as with the previous example we can see clear differences in the evolution of coefficients between the simplified Maxwell and Maxwell potential in Fig. 16.     [9] (reference, solid lines), with values kindly provided by the authors. Note, that we transformed our coordinates to match the orientation of the reference values for this plot When looking at the heat flux and stress tensor, we find good correspondence to previously published results in [9] in Fig.17, with an absolute difference in the order of 10 −6 . Note, that the approximation space for this comparison is not the same. Referring to the method demonstrated in [9], our ansatz corresponds to M 0 = M = 20, while the values provided by the authors have been calculated with M 0 = 15, M = 60. The good correspondence suggests that this difference has little influence on the moments of this low order and that the irreducible Burnett ansatz produces results very similar to the spectral Burnett ansatz. This should be the case, since the difference between the spectral and the irreducible Burnett ansatz lies solely in the exploitation of the collision tensor's decomposition.

Summary
With the irreducible Burnett ansatz we were able to reduce the computational effort of calculating Boltzmann's collision operator in the spectral Burnett ansatz both in computational time and memory by implementing an algorithm inspired by the representation theoretical decomposition of the collision operator. This decomposition relies on the fact that the used basis is adapted to the irreducible subspaces of the approximation space.
This decomposition not only allowed us to find an efficient algorithm, but also provides an understanding behind the structure of the collision tensor by separating factors that depend on the underlying physics from those that are of pure mathematical nature: The bilinear collision tensor can be expressed as product of the impact tensor with coefficients depending on the potential between the colliding particles and the coupling tensor with coefficients only depending on the basis choice of the irreducible subspaces.
With spherical harmonics as basis, further structure of the coupling tensor could be exploited by looking at O 2 -irreducible subspaces. For example, the coupling tensor could be split into the case where all anisotropic indices are zero and everything else. This is not necessarily generally applicable, but there may exist interesting use-cases.
Applying the method of this paper to collision operators of different kind, such as polyatomic gases [28] or the Landau operator [29], looks like a promising next step. not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
and similarly for r · c 1 : We further note, that g and n are affected by the variable change: Note how the transformation r acts on the normal vector n: The matrix R now rotates around the axis r (n), instead of the axis n, while the integration domain is transformed by r −1 . Therefore, we can equivalently integrate over n ⊥ r c −c 1 and substitute all n by r −1 n. Now we are left with an expression for S(r −1 ( f )) that has r (c) instead of c, andc 1 instead of c 1 everywhere, which is equivalent to the right hand side of (50). This shows covariance of the collision operator. Note the core reason for this: The collision operator is a concatenation of rotationally covariant operations such as taking the norm of a vector, or integrating over all directions.

Representation Theory: An Overview
Below we describe the core aspects of representation theory that we employed in this paper. This is a use-case restricted overview of the subject and we refer the reader to lecture books, e.g. [30][31][32][33][34][35], and to the original material such as [17,36] for further reading and full proofs of the following statements. Definition 9.1 (Representation (Definition 4.1 from [33])) Let G be a group and F a field. A representation of G over F consists of called the operation of the group on the vector space, that satisfies the axioms in Table 5.
We can say V is a representation of G and G operates on V .
For our purposes, we are interested in representations of the orthogonal groups, mainly O 3 given by orthogonal matrices r ∈ R 3×3 , i.e. G = O 3 . The field F is R and the vector spaces V we use are the tensor space of R 3 and the space of real polynomials in three arguments R[x, y, z]. For both the operation of O 3 is given in terms of the matrix vector-product r · v with v, v i ∈ R 3 : Non-triviality or normalization ∀v ∈ V : Axioms of homomorphisms Meaning One can generalize this to the tensor product between two vector spaces V , U each a rep- For two representations to be isomorphic, their vector spaces also need to be isomorphic. If the dimensions of two vector spaces are not equal, then they cannot be isomorphic representations.
Since the collision operator describes a physical process, one of its fundamental properties is, that it does not depend on the choice of coordinates. We also note, that it is a bilinear map between polynomial spaces, C : P M × P M → P M , which can be written as a linear map between the tensor product of two polynomial spaces to a polynomial space, P M ⊗ P M → P M . In mathematical terms, the collision operator must be a homomorphism from Hom RO 3 (P M ⊗ P M , P M ), i.e. a natural map. While this already restricts the degrees of freedom when compared to general linear map P M ⊗ P M → P M , we can go further by using the concept of irreducible subspaces. [33])) A subspace U of any vector space V is called invariant with respect to the operation of a group G, or G-invariant, iff for all g ∈ G and for all u ∈ U , the operation of g on u still is in U , g · u ∈ U or equivalently: ∀g ∈ G : g(U ) = U . Definition 9.4 (Irreducible subspaces (Definition 4.2 from [33])) A subspace U = 0 of any vector space V is called irreducible with respect to the operation of a group G, or Girreducible, iff it is invariant and there is no invariant subspace W ⊂ U with the properties W = U and W = 0. If a subspace U is invariant and dim U = 1, then it is also irreducible, however an irreducible subspace can be of dimension larger than 1.

Definition 9.3 (Invariant subspaces (Definition 4.2 from
As mentioned in [17], the space of symmetric trace-free tensors of order n, STF n (R 3 ), is O 3 -irreducible. With the isomorphism between STF n (R 3 ) and H n , which is also explored in [17], H n is also irreducible. Remark 9.5 Two irreducible subspaces U 1 , U 2 that are not isomorphic to each other are orthogonal to each other with respect to any scalar product of V that is compatible with the operation of the group G, i.e. is itself a natural map: Now we can decompose our bilinear collision operator as the sum over natural maps between irreducible subspaces, and using the isomorphism between H n a and H n , we can restrict the different types of maps to We now note that symmetric tensors are isomorphic to polynomials of the same degree and same underlying vectorspace, Chapter B.2 in [38], i.e. with R 3 we can map symmetric tensors of degree n from the space of symmetric tensors Sym n (R 3 ) to homogeneous polynomials of degree n in R[x, y, z]. We identify the tensor space that is isomorphic to the space of harmonic polynomials H n of degree n as the space of symmetric trace-free tensors STF n of same degree, as is explored in [17]. Anything that we learn from the tensorial view, we can apply to the polynomial view. We first note, that the tensor product of two representations is itself a representation. We also note, that the codomain of the natural map is an irreducible representation of O 3 , so using Schur's lemma, we know that the natural maps in (66) can only be non-zero if the codomain is isomorphic to one of the irreducible subspaces of the domain. Finding the irreducible subspaces of the tensor product of two representations is not trivial in general and this decomposition is called Clebsch-Gordan theory [39]. For the case of O 3 and symmetric trace-free tensors over R, one can find the decomposition of the maps in (66) relatively easily and we show the argument only for the sake of completeness. We start by considering the types of natural maps between tensors [36]: • Permuting tensor factors, • Taking the trace between two factors, • Taking the tensor product with the equivalent of c 2 , ω := i b i ⊗ b i for any orthonormal basis set { b 1 , b 2 , b 3 } of R 3 , also known as the δ i j -tensor.
Any natural map between tensor spaces can be expressed as linear combination and composition of the above maps. They can be represented graphically with Brauer diagrams [36]. This gives us a search space for the isomorphism STFn ⊗ STFñ → STF n . We first note, that the result must be a tracefree tensor, so we use the detracer D from [17]. D is an explicitly given projection from a symmetric tensor onto the symmetric tracefree tensor space of same order. As its input is a symmetric tensor, we need the symmetrizer (Part B.2 in [38]) which gives us the overall projector D • q, which is applied last. From Remark 9.5 and because the multiplicity of STF n in Sym n is one, we know that this map is unique up to a constant. What remains to be shown is for any natural linear map φ that we apply before this projector such that the overall map D • q • φ is STFn ⊗ STFñ → STF n , the overall maps at most differ by a constant factor and thereby fulfilling the promise in Remark 9.7. In other words: We want to show that span D • q • φ : STFn ⊗ STFñ → STF n φ is a natural linear map is one-dimensional or the space of maps to zero. Next, we note that the overall tensorial degree of the domain isn +ñ, and we would like to show the restriction n ≤n +ñ. To do so, we find all natural maps STF n → (R 3 ) ⊗(n+ñ) and show that they must all be zero, leaving us with the conclusion that a representation isomorphic to STF n cannot be contained in (R 3 ) ⊗(n+ñ) . Any natural map STF n → (R 3 ) ⊗(n+ñ) must reduce the tensorial degree n ton +ñ, leaving the trace as the only candidate, but since the tensors in STF n are trace-free by definition, the result is zero. First taking the tensor product, i.e. multiplying, with any number of ω also does not yield a different result, since in the end more traces need to be taken than the number of ω that have been multiplied, inadvertently leading to taking at least one trace of the tensor in STF n . Therefore, for n >n +ñ there does not exist a subspace in (R 3 ) ⊗(n+ñ) isomorphic to STF n , leading us to the desired restriction.
We further note, that the maps listed above all change the degree of tensor product by an even number, soñ +n − n ∈ 2N 0 . It is clear that the traces can only be taken between one factor of STFn and one factor of STFñ each, as opposed to a trace of two factors in STFn, which would map to zero by definition of STFn. This gives the third condition on the anisotropic indices from (32),ñ +n − n ≤ 2 min(ñ,n). Finally, we note that the codomain is also tracefree, so if somewhere in the composition the tensor product with ω is taken, an additional trace needs to be taken as well to achieve tracefreeness, otherwise the detracer would project onto zero due to the orthogonality between tracefree tensors and tensors with a trace (Part B, section 6 in [40]). Overall, the space of homomorphisms STFn ⊗ STFñ → STF n is one-dimensional, so each map is unique up to a constant.
Following the isomorphism back to the polynomial spaces, we can apply the results oneto-one, leading to the unique up to a constant natural map given by the coupling tensor. Looking at (64): and picking one non-zero instance for each of the maps that satisfy (32): Noting that due to the isomorphism between the irreducible subspaces, allQ nnñ abc with same combination of anisotropic indices work the same, we arrive at the decomposition of the collision tensor from Eq. (33).