Using Spherical Harmonics to solve the Boltzmann equation: an operator based approach

The transport of charged particles or photons in a scattering medium can be modelled with a Boltzmann equation. The mathematical treatment for scattering in such scenarios is often simplified if evaluated in a frame where the scattering centres are, on average, at rest. It is common therefore, to use a mixed coordinate system, wherein space and time are measured in a fixed inertial frame, while momenta are measured in a"co-moving"frame. To facilitate analytic and numerical solutions, the momentum dependency of the phase-space density may be expanded as a series of spherical harmonics, typically truncated at low order. A method for deriving the system of equations for the expansion coefficients of the spherical harmonics to arbitrary order is presented in the limit of isotropic, small-angle scattering. The method of derivation takes advantage of operators acting on the space of spherical harmonics. The matrix representations of these operators are employed to compute the system of equations. The computation of matrix representations is detailed and subsequently simplified with the aid of rotations of the coordinate system. The eigenvalues and eigenvectors of the matrix representations are investigated to prepare the application of standard numerical techniques, e.g. the finite volume method or the discontinuous Galerkin method, to solve the system.


INTRODUCTION
The transport of charged particles, photons or neutrinos in an inhomogeneous scattering medium is a recurring problem that appears across multiple sub-fields of physics (e.g.Braginskii 1958;Chandrasekhar 1960;Ginzburg & Syrovatskii 1964;Mihalas & Mihalas 1984;Lewis & Miller 1984).Astrophysical contexts include neutrino and radiation hydrodynamics in supernova explosions, as well as the propagation and acceleration of energetic charged particles (cosmic rays).Typically, one considers the evolution of the single particle phase-space density,  (x, p, ) by solving a transport equation which can be written in the form d  d = Ĉ  .Here  is an affine parameter that defines a point on a photon or particle trajectory (for a photon or neutrino the distance along a ray, for a particle with mass the proper time) and Ĉ is a collision operator that describes stochastic disturbances of these trajectories.For photons and neutrinos, the trajectories are simply geodesics in the local space-time, and Ĉ contains the scattering, absorption and emission cross-sections.These may be prescribed in the scattering-centres' rest frame, which may be determined by solving a separate set of equations (e.g., those of hydro-or magnetohydrodynamics).These equations in turn can be augmented by terms containing moments of  , i.e. feedback.In general, the rest frame of the scattering centres is non-inertial, and this must be taken into account when formulating the relevant transport equation (e.g.Achterberg & Norman 2018).
For example, the phase-function for photon scattering is usually specified in terms of the coordinates of the momenta p and p ′ of the photons involved, measured in the "co-moving" frame in which the scattering centres are locally at rest.This leads to a formulation of the radiation transport equation in so-called "mixed" coordinates (Lindquist 1966;Castor 1972;Riffert 1986), with the phase-space density treated as a function of position x and time , both measured in the laboratory frame, while the momentum p, is measured in the co-moving frame.
For cosmic rays, the situation is more complicated.As in the case of photons and neutrinos, one conventionally assumes the presence of a "background" plasma whose evolution can be determined, for example via the equations of ideal magnetohydrodynamics. But, in this case, these are augmented not only by terms in the stress-energy tensor, but also by the current and charge density carried by the cosmic rays.The charged particle trajectories are then not geodesics, but are determined by the electromagnetic fields carried by the background plasma.In dense plasmas, where Coulomb collisions among the different charged species occur, Ĉ can be modelled in the co-moving frame as a Fokker-Planck operator (Rosenbluth et al. 1957;Shkarofsky et al. 1966).However, in dilute, astrophysical plasmas, Coulomb collisions are unimportant, and Ĉ is instead used to model the effects of scattering on stochastic turbulent fluctuations that are superimposed on the large-scale fields carried by the background plasma.Given suitable constraints on the properties of this turbulence, Ĉ can again be formulated as a Fokker-Planck operator, this time in the (non-inertial) wave-frame (e.g. Skilling 1975;Webb 1985;Kirk et al. 1988) 1 .In principle, the large scale fields take into account the current and charge density of the cosmic rays.Consequently, the left-hand-side of the transport equation resembles a Vlasov equation, which has given rise to the name "Vlasov-Fokker-Planck equation".
For both photons and charged particles, the single particle phase-space density is a function of six independent variables and time, making direct discretisation approaches prohibitive.In the case of the space and time variables, numerical techniques for solving the MHD equations are well-known.These can be applied also to the VFP equation, since  is not expected to vary on scales shorter or faster than those of the MHD variables.Also, CR distributions are generally smooth functions of the particle energy (or the modulus  of the momentum), so that discretization of this variable is not problematic.However, the two remaining angular variables (in momentum space) pose a problem, since strong anisotropies can develop, in particular close to oblique shock fronts in the background plasma or in the presence of highly relativistic flows.An effective approach here is to expand the phase-space density in a series of known functions, and thereby resolve the anisotropy by extending the expansion to higher order until convergence is achieved (Kirk & Schneider 1987;Kirk et al. 2023).In the case of cosmic rays in non-relativistic flows, a suitable set of functions is that of spherical harmonics.This is the approach we adopt here, by expanding the single particle distribution in a series of spherical harmonics    , which are functions of the spherical polar coordinates  and  of the momentum vector p (defined with axis along   ): or, equivalently, using Cartesian tensors (Schween & Reville 2022).The equivalence is particularly well exploited in the (projected-)symmetric trace-free formalism described by Thorne (1980Thorne ( , 1981)).Alternatively, it is possible to derive the equations for the moments of the distribution function and, thus, to directly investigate the relevant transport coefficients, see for example Webb (1989) and references therein.Expanding the single particle distribution function reduces the number of dependent variables from six to four, though it requires a procedure to determine an effective transport equation for each coefficient in the expansion.While derivations of this system of equations can be found in the literature for different applications, both for Cartesian expansions (see Johnston 1960;Williams & Jokipii 1991;Thomas et al. 2012, truncated at low order) and spherical harmonics (e.g.Bell et al. 2006;Tzoufras et al. 2011;Reville & Bell 2013), the resulting equations are cumbersome, and the physical meaning is obfuscated.
In this paper we present a simplification of this expansion procedure which avoids lengthy algebraic manipulations, when expanding to high order.This facilitates checking the convergence of the expansion and allows for a direct physical interpretation at each stage.We restrict our derivation to the propagation and acceleration of charged particles in a prescribed background plasma.In Section 2 we introduce the relevant transport equation, namely a variant of the Vlasov-Fokker-Planck (VFP) equation, and simplify it by making the additional, although not necessary, assumption that the co-moving frame moves slowly in the laboratory frame, and that the scattering operator is isotropic.For such a scattering operator, the eigenfunctions are the spherical harmonics used in the expansion.This motivates an approach which exploits operators acting on the space of spherical harmonics and their matrix representations.In Section 3 we identify and define the operators which are present in the VFP equation and show that the system of equations can then be obtained by replacing the operators with their matrix representations.This has the advantage of retaining the original structure of the VFP equation.We explicitly compute the matrix elements of the matrix representations.Moreover, in Section 4 we prove that the effort of computing them can be considerably reduced with the help of rotations of the coordinate system in momentum space.We investigate how these rotations act in the space of spherical harmonics and how they change the operators and their matrix representations.The key insight is that the operators' matrix representations can be considered as "rotated" versions of each other and thus it is possible to compute any one of them by rotating another.We quote explicit expressions for the necessary rotation matrices.Having derived and simplified the computation of the system of equations, we note that it contains redundant equations.The reason is that the reality of the single particle phase-space density imposes symmetries on the expansion coefficients.In Section 5 we remedy this using the real spherical harmonics.The relation between the complex and real spherical harmonics can be expressed with the help of a basis transformation matrix.Multiplying the system of equations with this basis transformation matrix removes the superfluous equations.Numerically solving the system of equations using standard techniques, e.g. the finite volume method or the discontinuous Galerkin method, requires knowledge about the eigenvalues and eigenvectors of the matrix representations.In the last part of this paper, i.e . in Section 6, we prove that the representation matrices have the same eigenvalues and that there eigenvectors are also rotated versions of each other.Moreover, we show that eigenvalues are the roots of the associated Legendre polynomials.The implementation of the discontinuous Galerkin method exploiting this knowledge is the subject of a companion paper.

PARTICLE TRANSPORT IN PLASMAS
The transport of charged particles and photons in a scattering, non-stationary medium was recently revisited by Achterberg & Norman (2018).Their work provides a rigorous derivation of the mixed-frame transport equation in the context of special relativity.We adopt the mixed-frame transport equation as our starting point, and refer the reader to their paper for further details.
We consider particles that move in a background fluid, whose velocity field can be described with the 4-vector   = Γ(, U) (bold symbols denoting 3-vectors), with Γ = (1 −  2 / 2 ) −1/2 the bulk Lorentz factor of the background fluid.For the particles on the other hand, the four-velocity is   = (, V), with  = / = / 2 the associated particle Lorentz factor.The background plasma is permeated by electric and magnetic fields, E and B. Using unprimed quantities to denote quantities measured in the laboratory frame, and primed quantities to denote those measured in the fluid rest frame, one arrives at the following transport equation where, on account of the mixed coordinates, the proper-time derivative along the particle trajectory is (2.2) As described in Achterberg & Norman (2018, eq. 15), the second two terms inside the square brackets of equation (2.1) correspond to the fictitious forces introduced from the accelerating frame.The ω ℓ vector stems from changes in the direction of U, with l a unit vector along U =  l .
The particular form of the collision term we adopt here corresponds to isotropic scattering in the local fluid frame with rate  ′ (  ′ ), although such a simplification is not essential.The symbol Δ  ′ ,  ′ denotes the Laplace operator in spherical angles, (2.3) The collision term thus assumes that the scattering is elastic, changing only the direction of a particle's momentum and not its magnitude which need not be the case if the scattering centres are dispersive, or anisotropic (Melrose 1969;Skilling 1975).A detailed study of these effects is beyond the scope of the current paper.
For  ≪ , we can expand equation (2.1) retaining all terms to first order in /: (2.4) In the rest of this text we will work with eq.(2.4), but we emphasise that the method to derive the equations for the expansion coefficients of the spherical harmonics can be straightforwardly applied to the fully relativistic VFP equation (2.1).We will henceforth drop the primes from all co-moving terms, and it is to be understood that U, V, p, , B and  are referring to quantities defined in the rest frame of the thermal plasma.We also restrict our attention to ideal MHD, for which E ′ = 0, though since it independent of V, it can be re-introduced if required by adding to the U/ term.

DERIVATION OF THE SYSTEM OF EQUATIONS
We now convert eq.(2.4) into a matrix equation by inserting the expansion (1.1) and projecting onto the spherical harmonics.Projection involves integrals over the unit sphere.We adopt a shorthand notation for the corresponding scalar product: This yields (3.2) Eq. (3.2) can be interpreted as a matrix equation.The scalar products are the elements of matrices that multiply the vector ( f )  (,)    or its derivatives.The index  used here is a one-to-one function of the indices  and .We call this function an index map and it determines how the expansion coefficients are ordered.An example order is given through f = (  0 0 ,  1 1 ,  0 1 ,  −1 1 , . . . )  .In this case the index map  is  (, ) =  ( + 1) − , with  starting at zero.Its (unique) inverse is given as  = ⌊ √ ⌋ and  =  ( + 1) − .

The identity operator and the collision operator
As a first example, we consider the time derivative.This can be written in the form: where 1 is the identity matrix.Note that we introduced a second index map ( ′ ,  ′ ) =  ′ ( ′ + 1) −  ′ such that the notation resembles ordinary matrix-vector product.In particular, we emphasise that summing over  and  can be interpreted as the multiplication of the identity matrix with the time derivative of the vector f .Each time additional matrix indices are required a "copy" of the index map  (or ) will be used, see for example the matrix-matrix products in eq.(3.42).
The identity matrix can be considered to be the matrix representation of the identity operator 1 acting on the space of the spherical harmonics.We define the matrix representation of an operator Ô acting on an inner product space H = span{ 1 ,  2 , . . .} to be where the   s are the basis vectors of H . (3.4) We note that the space of all spherical harmonics is a Hilbert space, whose inner product is defined by the orthogonality relation We denote this space by The collision operator serves as another example for the interpretation of the integrals as matrix elements.The right-hand side of eq.(3.2) is Where C is the matrix representation of the collision operator in the spherical harmonic basis.It is a diagonal matrix and its elements are easily computed, since the spherical harmonics are eigenfunctions of the angular part of the Laplace operator, see eq. (2.3).They are These two examples suggest that the system of equations which determine the expansion coefficients    can be formulated as representation matrices of operators acting on the space of spherical harmonics.We thus focus on identifying other operators whose action on the spherical harmonics are either known or easily derived.

The angular momentum operator
The magnetic force term in the VFP equation (2.4) is the counterpart in momentum space of the angular momentum operator in configuration space that is well-known from quantum mechanics: L QM −ix × ∇  .I. e., where ω = B/ is the angular frequency vector, and L = −ip × ∇  .We choose to keep the name "angular momentum operator" for L, despite its action on the momentum space variables p In terms of the coordinates  and , The operator L± is the angular momentum ladder operator.These definitions can for example be found in Landau & Lifshitz (1977, eq. 26.14-15).We express L and L in terms of these ladder operators and we get for the angular momentum operator We stress that the angular momentum operator acts only on  and .Hence, the magnetic force term in eq.(3.2) becomes where we took    out of the integral and introduced the representation matrices   of the angular momentum operators.Note that we implicitly sum over the index  ∈ {, , }.We will frequently write   to refer to all three representation matrices at once.
The action of the angular momentum operators on the spherical harmonics is well-known from quantum mechanics.For example, in Landau & Lifshitz (1977, §27) their action is given to be (3.13) The matrix representations of the angular momentum operators can now be computed using eq.(3.13) and eq.(3.11).They are Alternatively, these matrix elements can also be found in Varshalovich et al. (1988, Chapter 13.2, eq. 42).

The direction operators
In addition to the identity and angular momentum operators, another set of operators appears in the VFP equation.We call them direction operators because they "point" in the direction of the coordinate axes: (3.17) These operators act on the space of spherical harmonics, i.e.Â : S → S. They appear, for example, in the part of the VFP equation we call the spatial advection term.In the equations for the expansion coefficients (3.2), the spatial advection term is (3.18) They also appear in front of the time derivative of the distribution function, i.e. 1 As before, the computation of the matrix representations of the Â operators leads us to investigate their action on the spherical harmonics.There are multiple ways to determine their action.For example, it is possible to use recurrence relations for the associated Legendre polynomials or to use ladder operators, which increase/decrease  and/or , to represent them.
We begin with Â and we use a recurrence relation to examine its action, namely With its help, we show that the action of Â is (3.21) Instead of using a recurrence relation for the associated Legendre polynomials, we can work with ladder operators which increase (decrease) .These operators are investigated in Fakhri (2016, eq.12) and they are defined as Ĵ± () ± sin  +  cos  .(3.22)With this definition at hand, a short computation shows that The action of Â is now given by the action of the ladder operators Ĵ± , which can as well be found in Fakhri (2016, eq.14 a,b), and it is For the sake of completeness, we proceed analogously with the operator Â and Â , i.e. we, first, compute their action using a suitable recurrence relation for the associated Legendre polynomials and, secondly, we express their action with the help of appropriate ladder operators.But already at this point, we highlight that an explicit computation of the representation matrices of Â and Â is not necessary, as they are connected via rotations (see section 4.2).
To compute the action of Â and Â , we use the two recurrence relations (3.26) With these relations at hand, we show that the action of Â is The coefficient    is In a completely analogous manner the action of Â can be shown to be Again, the action of Â and Â can alternatively be derived with ladder operators which can be found in Fakhri (2016, eq. 25, eq. 34), i.e.
The square brackets denote the commutator of two operators, i.e. [ Â, B] Â B − B Â.These operators shift  and  either in the same "direction" or in opposite "directions".Their action is described by Fakhri (2016, eq.26 a,b, eq.35 a,b) and it is (3.32) The operators Â and Â can be represented with the help of the above ladder operators, i.e.
This can be verified by direct computation.Knowing the action of the direction operators (see eq. (3.21), (3.27) and (3.29)) , we can now compute their matrix representations.This yields These matrix elements can also be found in Varshalovich et al. (1988, Chapter 13.2, eq.14-16).

Products of operators
The only term in the VFP equation that cannot immediately be formulated using the identity, angular momentum and direction operators is This term is a consequence of the transformation of p to the local fluid frame and is referred to as a fictitious force term (Achterberg & Norman 2018).
We show now that ∇  can be expressed in terms of the direction operators Â and the angular momentum operators L : Here, e  , e  and e  , are unit vectors of the spherical coordinate system, and we point out that the components of e  are just the direction operators: Note that expressing ∇  with the help of the direction operators and the angular momentum operators comes at the cost of introducing products of operators.Up to now, we saw that it was possible to derive the system of equations for the expansion coefficients by simply replacing the operators with their respective matrix representations.We now prove that this is still true for products of operators if some care is taken.
First of all, we show that the matrix representation of a product of operators is the product of the matrix representations of the operators.To see this, we point out that the identity operator in the space of the spherical harmonics can be written as If  ∈ S, then it can be written as a linear combination of spherical harmonics and the scalar products in the above sum yield the coefficients of this linear combination. 2or the matrix representation of the product of two operators, say Â and B, we get (3.42) Hence, the matrix representations of the two operators are multiplied.We emphasise that the above statement does not hold if we truncate the expansion of the distribution function at a finite  max .In this case, we compute the matrix representations of the operators acting on the finite dimensional space S  max span{ 0 0 ,  is not an element of the space S  max , it makes sense to define that the result is zero when they act on    max .But if a product of operators appears whose first operator increases , while the second operator decreases it again, then the above definition would imply that this product yields zero when it acts on    max .This is not the correct result, because the joint action results in a constant times    max , which can be represented in S  max .This difficulty can be avoided, if we compute the matrix representations of the operators for  =  max + 1, evaluate the products and, subsequently, reduce the matrices to  max .
Since we now know how to compute the matrix representations of a product of operators, we can now look at eq. (3.2) and see how the fictitious force term contributes to the system of equations, i.e. (3.43) We stress that it is exactly this computation where the proposed operator method provides the most obvious benefit.The repeated application of recurrence relations for the associated Legendre polynomials is replaced with the evaluation of matrix products, which avoids lengthy and cumbersome calculations (see Reville & Bell 2013, Appendix A).

The complete system of equations
We conclude this section with a presentation of the complete system of equations, which determines the expansion coefficients    , and an explanation of its terms.The system of equations is (3.44) The first term is the time derivative of the expansion coefficients including a correction of time stemming from the transformation of momentum variables into the rest frame of the magnetised (or background) plasma.Its second term models how the expansion coefficients change because of the motion of the plasma, which is described by U, and the motion of the particles given through A  .As already stated above, we refer to this term as the spatial advection term.The next term changes the energy contained in the expansion coefficients because of the fictitious force acting on the particles.We call it the momentum advection term, because it advects the coefficients in -direction.The fourth term reflects that the fictitious force not only accelerates (or decelerates), but it also changes how the energy is distributed among the expansion coefficients.The fifth term represents the effect of the magnetic force, namely the rotation of the velocity of the particles.This rotation translates into a rotation of the expansion coefficients. 3We refer to this term as the rotational term.The last term is modelling the interactions of the particles with the thermal plasma.C is a diagonal matrix, which causes the expansion coefficients with  ≥ 1 to decay exponentially at the rate  ( + 1)/2.The coefficients with  ≥ 1 model the anisotropies of the distribution function and the effect of the exponential decay is to drive the distribution towards isotropy.We call this term collision term.
We, once more, would like to highlight the one-to-one correspondence between operators and representation matrices.The system of equations can be derived by simply replacing the operators with their corresponding matrices.

ROTATIONS IN THE SPHERICAL HARMONIC SPACE
In this section we present an alternative way to compute the representation matrices A  , A  ,   and   .It is based on the observation that our choice of the coordinate system in momentum space was arbitrary and that the direction operators Â and the angular momentum operators L are defined with respect to the coordinate axes, see their definitions in eq.(3.17) and eq.(3.10) respectively.The spherical harmonics are also defined in the same coordinate system.Rotating them clockwise by /2 about the -axis is equivalent to defining them in an equally rotated coordinate system.In the rotated system the former -axis is the -axis and we expect the representation matrix A  to have transformed to A  .Taking advantage of this requires us to investigate how to rotate the spherical harmonics and how the operators and their representation matrices change under such a rotation.

Deriving the rotation operator
We begin our investigation with formally introducing a rotation operator in the space of spherical harmonics, i.e.  m  (, ) = R   (, ).m reflects the idea that we can think of the rotated spherical harmonics as being defined with respect to a new polar direction given by a corresponding rotation of the original one.A direct consequence of this thought is that a rotation does not change the degree  of the spherical harmonic.
An explicit form of the rotation operator R can be derived by defining its action on the spherical harmonics and by studying an infinitesimal rotation.This is, in essence, an application of representation theory.We present the details here for completeness.In doing so, we follow Jeevanjee (2011), in particular Chapter 5 and example 5.9 therein.
We point out that and define the action of R to be Where R is a rotation matrix which rotates the vector p.
We now look at an infinitesimal rotation to derive an explicit expression for the rotation operator.To simplify the discussion, we restrict ourselves to a rotation about the -axis by an infinitesimally small angle d, i.e.R  (d)p = p + de  × p.Its inverse is R −1  (d) = p − de  × p.An explicit expression for the rotation operator is now constructed by plugging the infinitesimal rotation into the definition of its action (4.1).This yields Where we Taylor expanded in d at zero.In a next step, we set d = / and perform  infinitesimal rotations, which add up to a rotation about the -axis by an angle .In the limit  → ∞, we get We state, without further proof, that this generalises to R =  −ir•L for an arbitrary rotation about the unit vector r by an angle .
The matrix elements of the rotation operator R are the Wigner-D functions, which are well-known in the quantum theory of angular momentum, see Varshalovich et al. (1988, Chapter 4.5, eq. 1).They are defined to be (4.4) The conjugate transpose of the Wigner-D matrix is where we used that for arbitrary matrices X and Y the following two properties hold: Firstly, for (X • • • X) † = X † • • • X † the order of the matrices does not matter and, second, that (X + Y) † = X † + Y † .Moreover, we took advantage of the fact that the representation matrices of the angular momentum operator are Hermitian, namely   † =   , see eq. (6.3, 6.4).
A direct and useful consequence is that the Wigner-D matrix is unitary4 , i.e.

The change of operators and their representation matrices under a rotation of the coordinate system
Knowing how to rotate the spherical harmonics, we are in a position to compute the action of an operator Ô in the rotated coordinate frame.
We once more note that we can compute the spherical harmonics defined in the rotated frame simply by rotating the original spherical harmonics in the same manner as the coordinate system.
Before we derive the action of an operator in the rotated coordinate system, we would like to make a technical remark: Transformations of coordinates can be interpreted actively or passively.An active rotation, for example, yields the coordinates of a rotated vector, whereas a passive rotation gives the coordinates of the vector in a rotated coordinate frame.Active and passive rotations are equivalent if the sign of the rotation angle is changed, i.e.R pas (r) = R act (−r).Even though we write as if we rotated the coordinate system, we actually interpret the rotations as active.Instead of working with two different coordinate systems in which the spherical harmonics are defined, we work with two different bases of the space of spherical harmonics, namely the original spherical harmonics and the rotated spherical harmonics.Both are defined in the original coordinate system.The reason to take the active point of view is to allow for a direct application of representation theory without the difficulties from two different sets of coordinates, for example, p in the rotated and p in the original coordinate frame.
That said, let Õ be the operator in the rotated coordinate system, then Õ m  (, ) = R Ô R−1  m  (, ); the inverse rotation R−1 "brings" the spherical harmonic back into the unrotated coordinate system.The action of Ô on the unrotated spherical harmonic is known.Its result is rotated again and Õ gives, as expected, results in the rotated coordinate system.
Eventually, we are interested in the representation matrices of operators and, thus, we would like to know how the representation matrix of an operator Ô changes under a rotation of the coordinate system.The representation matrix Õ in the rotated frame is given by the action of Ô on the rotated spherical harmonics, i.e. (4.7) Or compactly, Õ = U † OU.

Computing the representation matrices using rotation matrices
As explained in the introduction of this section, the operator Â has to be Â in a coordinate system which is rotated by −/2 about the -axis.
If we know the representation matrix A  , we can compute how it looks like in the rotated frame using the formula in eq.(4.7) and we know that in this frame it has to equal A  .Hence, Moreover, we compute the representation matrix A  knowing A  by rotating the original coordinate system about the -axis by −/2, i.e.
This leads us to the conclusion, that it is enough to compute the representation matrix of Â .We highlight that we do not have to compute the representation matrices of the involved rotations ourselves.Since their matrix elements are the Wigner-D functions, they can be, for example, be found in Varshalovich et al. (1988, Section 4.5.6 eq. 29, eq. 30 and Section 4.3.1, eq. 2).They are where the limit of the sum is  = min( −  ′ ,  − ).
We remark that the representation matrices of the angular momentum operators can as well be obtained by rotating   , i.e.
We emphasise that to construct the system of equations (3.44), we only need four matrices, namely A  ,   , U(−/2e  ), and U(−/2e  ) and two of them are diagonal matrices.This should considerably diminish the burden of its implementation in future numerical applications.

A REAL SYSTEM OF EQUATIONS
The phase-space density  is a real function, i.e  ∈ R.This implies the existence of a relation between the coefficients    of the spherical harmonic expansion of  as given in eq.(1.1), because the results of the sums must yield a real function.Since a scalar (function)  is real if and only if  * =  , the relation between the coefficients can be shown to be where we used that    * = (−1)   −  .This means that the system of equations (3.44) contains too much information, i.e. too many equations, because it determines the expansion coefficients  −  as well.And these, as we have just proven, can be directly obtained from the ones with positive .A way to remedy this is to use real spherical harmonics instead of complex spherical harmonics.We note that the above derivations are greatly simplified by working with the complex spherical harmonics, and this is why we choose to transform to the real spherical harmonics only after the matrix representations have been determined.

Definition of the real spherical harmonics and the basis transformation
We begin with a definition of the real spherical harmonics.One way to arrive at their definition is to rewrite the spherical harmonic expansion (1.1) as (5.2) In the last line we defined the real spherical harmonics   .The expansion coefficients   are real as well.Furthermore,  01 = 0 and  01 = 0.The real spherical harmonics are orthogonal functions too, i.e.
, which can be verified through a direct computation of the scalar product.
We can interpret the rewrite of eq. ( 1.1) as a basis transformation in the space of spherical harmonics S. The corresponding basis transformation matrix5 can be computed with the help of (5.3) We introduced a new index map  ( ′ ,  ′ ,  ′ ).It again encodes the ordering of the real spherical harmonics   and their coefficients   .For example, if the real spherical harmonics are ordered as S = span{ 000 ,  110 ,  100 ,  111 ,  220 ,  210 ,  200 ,  211 ,  221 , . . .}, the index map is  ( ′ ,  ′ ,  ′ ) =  ′ ( ′ + 1) + (−1)  ′ +1  ′6 .Note the index  starts at zero.The matrix elements of the basis transformation matrix S can be computed to be where we used the definition of the real spherical harmonics given in eq. ( 5.2) and, once more, the relation    * = (−1)   −  .Moreover, the basis transformation matrix has the useful property that it is a unitary matrix.This is a direct consequence of the fact that basis transformation matrices relating two orthogonal bases are unitary.We include a proof of this general statement for the case at hand, i.e. (1) (5.5) In matrix form this equation reads 1 = S † S and the inverse basis transformation matrix is S −1 = S † .

Turning the complex system of equations into a real system of equations
With the basis transformation and its inverse at hand, we can turn the complex system of equations (3.44) into a real system thereby removing the superfluous equations.To this end, we compute how the matrix representation of an operator changes under the basis transformation.We find for an arbitrary operator, say Ô, that Where the subscript  means that O  is the matrix representation with respect to the real spherical harmonic basis.Eq. (5.6) in matrix form reads O  = S † OS.
We now transform the complex system of equations (3.44) by multiplying it with S † from the left and inserting the identity matrix 1 = S † S where necessary.This yields (5.7) Note in particular that    = iS †   S includes the factor i. It is left to show that eq.(5.7) is actually a real system of equations, i.e. that all the appearing matrices and f  are real.We begin with the latter and show that the components of f  are (5.8) These are the expansion coefficients of the real spherical harmonic expansion, which we derived in eq. ( 5.2) and which, as we pointed out there, are real.To prove that the matrices are real, we had to compute them explicitly.The explicit expressions for all real representation matrices can be found in Appendix B.
At the end of this section, we stress that the real representation matrices can as well be computed with the help of rotation matrices, i.e. (5.9) (5.10)

EIGENVECTORS AND EIGENVALUES OF THE REPRESENTATION MATRICES
In the last section of this paper we investigate the eigenvectors and eigenvalues of the representation matrices of the direction operators and the angular momentum operators.The former play an important role in the implementation of numerical methods used to solve the system of equations, which we have derived.Examples are the finite volume method or the discontinuous Galerkin method.The statements about the eigenvalues of the representation matrices of the angular momentum operators are included for the sake of completeness.

General statements
A first important observation is that the matrix representations of the Â and L are Hermitian matrices, because the eigenvalues of Hermitian matrices must be real.Assume ξ is an eigenvector of an arbitrary Hermitian matrix A with eigenvalue  and with norm  √︁ ξ * • ξ = 1, then We now prove the above observation for the direction operators Â .We prove it for Â and only note that proof for the other two operators works analogously.We use the definition of the representation matrices given in eq.(3.18) to get the transpose conjugate representation matrix of Â , namely Hence A  is a Hermitian matrix.The same is true for the representation matrices of Â and Â .Hence, A  † = A  for  ∈ {, , }.
We proceed by proving that the representation matrices of the angular momentum operators are Hermitian.Knowing their action on the spherical harmonics, see eq. (3.13), we conclude that where we used that we get the same result when the roles of  ′ and  are interchanged.For L , we find (6.4) The proof for L is analogous.
A second statement about the eigenvalues of the representation matrices follows from the fact that the representation matrices are "rotated" versions of each other.This implies that their eigenvalues are the same.
The proof requires the concept of similar matrices: We say the matrices X and Y are similar, if they are related by Y = T −1 XT.Similar matrices have the same eigenvalues, because their characteristic polynomial is the same, i.e.

𝑝(𝜆)
det Eq. (4.9) says that A  is similar to A  and eq.(4.8) states that A  is similar to A  .The same is true for the representation matrices of the angular momentum operator, see eq. (4.11).Thus, all three matrices share the same eigenvalues.We define the spectrum to be the set of all eigenvalues of a matrix, namely (X) { ∈ C | ∃w ≠ 0 : Xw = w}.The fact that all three matrices have the same eigenvalues can now be expressed as (6.6) A third statement concerns the eigenvectors of the representation matrices: The eigenvectors are rotated versions of each other, i.e.
Where  is a diagonal matrix with the eigenvalues of A  on its diagonal and W are matrices whose columns are the corresponding eigenvectors.
An analogue computation can be performed for A  .

Eigenvalues of the representation matrices of the direction operators
In the previous section we showed that if we know the eigenvalues of one of the representation matrices, we know them for all.Since A  is simpler, i.e. it has less non-zero elements, we focus on determining its eigenvalues and, as we will see, its eigenvalues are the roots of the associated Legendre polynomials.
To see a general pattern emerge, we begin with computing the eigenvalues of A  for  max = 1 and  max = 2. Using the following ordering of the spherical harmonics for  max = 2, the matrix A  is, as we will see now, tridiagonal.We use eq.(3.35) to compute its matrix elements and introduce the shorthand    = [( + ) ( − )/(2 + 1)(2 − 1)] 1/2 for the coefficients in the referred equation.This yields for  max = 1 and for  max = 2 .We emphasise that the size of the block matrices A  , varies with  and that the pattern in eq.(6.8) is the same for arbitrary  max .The block diagonal structure of A  implies that its characteristic polynomial factors into the characteristic polynomials of its blocks, i.e.

𝑝(𝜆)
(6.10) Where we used that A  max , = A  max ,− , because    =  −  .A direct consequence of eq.(6.10) is that the eigenvalues of A  (i.e. the roots of its characteristic polynomial) are the eigenvalues of the block matrices    max, . We now prove that their eigenvalues are the roots of the th derivative of the Legendre polynomial  +1 .
Before we begin the proof we introduce (6.11) to denote the th derivative of the Legendre polynomial.Furthermore, we point out that    is part of the definition of the associated Legendre polynomials, i.e.    () = (−1)  (1 −  2 ) /2    ().Thus, if the eigenvalues of A  , are the roots of the th derivative of the Legendre polynomial  +1 , they are as well the roots of the associated Legendre polynomial   +1 modulo ±1.In a first step, we compute det(A  , − 1) developing it after the last row, i.e. (6.12) Eq. ( 6.12) can be read as a recurrence relation for the polynomials    .Fixing  to some integer value and setting   || () 1 the characteristic polynomials for arbitrary ,  can be computed recursively.
In a second step, we compare the recurrence relation (6.12) with the recurrence relation for the associated Legendre polynomials in eq.(3.20), which implies a recurrence relation for the th derivative of the Legendre polynomials, i.e.
−1 () .(6.13) We see that the recurrence relations in eq. ( 6.12) and eq.( 6.13) are similar.They differ in the factors in front of the polynomials on the right-hand side.This suggests that    () ∝    ().If the polynomials are proportional to each other, they have the same roots; multiplying a function with a constant does not change its roots.The factor of proportionality can be derived by noting that the leading order term of    () is (−1) −  − , using an explicit formula for the associated Legendre polynomials and dividing its leading order term by the factor which accompanies it.This results in We summarise and conclude that the eigenvalues of the representation matrices of the operators Â for a fixed  max are the roots of the polynomials    max +1 with 0 ≤  ≤  max , which are, modulo a constant of proportionality, the th derivatives of the Legendre Polynomial  +1 .Hence, the eigenvalues of A  are the roots of the associated Legendre Polynomials  0  max +1 , . . .,   max  max +1 modulo ±1.We note that the roots of the polynomials with  ≠ 0 have algebraic multiplicity two, see eq. (6.10).
The above conclusions have implications for the possible set of eigenvalues: First, the Legendre polynomials are orthogonal polynomials and the roots of orthogonal polynomials are real, simple and located in the interior of the interval of orthogonality, which is [−1, 1] in the case of the Legendre polynomials, see Milton & Stegun (1964, Section 22.16).Secondly, the roots of    with  ≠ 0 are the roots of the derivatives of the Legendre polynomials, i.e. the roots of    .Rolle's theorem states that the roots of the derivative of a continuous and differentiable function, must lie between the roots of this function.This implies that the roots of the polynomials    are "moving" towards zero with increasing , hence that all eigenvalues of A  are contained in the interval [−1, 1] and that the largest eigenvalue of A  is the largest root of the Legendre polynomial  0  max +1 =   max +1 .We can get a very good estimate for the largest eigenvalue using a formula to compute estimates of the roots of the Legendre polynomial of degree , which can be found in Tricomi (1950, eq.13).The formula is Note that the Legendre polynomials are even (odd) for even (odd) .The estimate for the largest eigenvalue of A  is  max ≈  1 .

Eigenvalues and eigenvectors of the sum of representation matrices
In the last part of this section we show how to compute the eigenvalues and eigenvectors of using the eigenvalues and eigenvectors of A  .In a first step we note that the sum of the representation matrices of the direction operators can be traced back to the scalar product η • e  =   Â +   Â +   Â , see eq. (3.40).Secondly, we define the unit vector n η/.Thirdly, we rotate the coordinate system such that the -axis of the rotated coordinate system is parallel to n.We denote the corresponding rotation matrix with R n .This changes the scalar product η (6.17)where we used that the coordinates ñ of n in the rotated coordinate system are e  .Eq. (6.17) shows that the sum of the direction operators given by the scalar product η • e  reduces to  Ã in the rotated coordinate system and so does its matrix representation (6.16) with respect to the spherical harmonics defined in it.
As explained in Section 4 transforming representation matrices into a rotated coordinate system requires to know how to rotate spherical harmonics.In eq. ( 4.4) we presented the necessary rotation matrix U(r) =  −i    , where the unit vector r is the axis of rotation and  is the angle by which the spherical harmonics are rotated.Because of eq.(4.7) and since we know that the sum of matrices (6.16) must equal A  in the rotated frame, The axis of rotation can be computed with r = e  × n/∥e  × n∥ and the angle  = e  • n.
We can now investigate the eigenvalues of the sum of the representation matrices (6.16).They are determined by whence we conclude that the eigenvalues of the sum of the representation matrices (6.16) are the eigenvalues of A  times , i.e. its th eigenvalue is where   is the th eigenvalue of A  .(6.20) We highlight that this relation between the eigenvalues of A  and the eigenvalues of the sum of the representation matrices is independent of  max and holds for each eigenvalue    .The geometrical interpretation of eq.(6.20) is that for a given value   the eigenvalues    lie on one of the sheets of a circular and two-sheeted hyperboloid which is parameterised by   and   .To illustrate this we plotted one eigenvalue of the sum of the representation matrices (6.16) for two values of   , see Fig. 1.
The eigenvectors of the sum of the representation matrices (6.16) are the rotated eigenvectors of A  , because We conclude this section with referring the reader to Garrett & Hauck (2016), who also investigated the eigenstructure of the A  matrices in the context of radiation transport and found similar results.Our work differs from theirs mainly in identifying the unitary matrices U as rotation matrices.Furthermore we would like to remark, that the eigenvalues of the real matrices A   are the same as the eigenvalues of the A  matrices, because S † A  S = A   is a similarity transformation.

CONCLUSIONS
The transport of particles or photons is very often modelled using variants of the Boltzmann equation.In many cases the scattering (or collision) operator has a particular simple form in the centre-of-momentum frame.This motivates to formulate the transport equation in a mixed-coordinate frame, i.e. a frame in which time and the configuration space variable x are measured in an inertial frame whereas p is given with respect to a non-inertial frame moving with the centre-of-momentum of the scattering centres.To render the equation amenable to be solved numerically, the single particle phase-space density is expanded in terms of known functions, for example the Legendre polynomials, to reduce the number of independent variables.This expansion leads to a system of partial differential equations evolving the expansion coefficients.
In this paper, a new method to derive this system of equations for the case of a spherical harmonic expansion of the single particle phase-space density has been presented.The method is based on operators acting on the space of spherical harmonics and their matrix representations.It was shown that the system is readily obtained by replacing the introduced and identified operators appearing in the VFP equation with their matrix representations.The final result is a compact matrix equation for the (real) expansion coefficients, which has the same structure and physical interpretation as the original VFP equation, i.e.
We simplified the computation of the matrix representations considerable with the help of rotations of the coordinate system in momentum space, see in particular eqs.(4.8), (4.9) and (4.11).We showed that the system of equations can be constructed from four representation matrices of whom two are diagonal and well-known.
The fact that the spherical harmonics (and the expansion coefficients) are complex whereas the single particle phase-space density is real, implies that the constructed system of equations contains linear dependent equations.We removed the redundant equations by multiplying the system with a basis transformation matrix, which turns the complex spherical harmonics into the real spherical harmonics, see eq. (5.7).
Eventually, we prepared the application of numerical methods, like the discontinuous Galerkin method or the finite volume method, through an investigation of the eigenvalues and eigenvectors of the representation matrices.A key insight consists in the fact that all representation matrices of the direction operators have the same eigenvalues.This holds true for the representations matrices of the angular momentum operators as well.Moreover, all eigenvalues are real.This was shown in Section (6.1).We emphasise, that we have proven that the eigenvalues of the representation matrices of the direction operators are the roots of the associated Legendre polynomials, see Section 6.2.We were also able to directly compute the eigenvalues and eigenvectors of the sum of the representation matrices of the direction operators, see eq. ( 6.20) and eq.(6.21).These results have all been exploited in the implementation of the Sapphire++ code, a discontinuous Galerkin scheme that solves equation (7.1).The details are described in a companion paper.
We conclude with two remarks: First, the presented method can be applied whenever a spherical harmonic expansion of the single particle phase-space density function is used.In particular, it is not necessary to formulate the transport equation in a mixed-coordinate frame.The only prerequisite is the ability to compute the representation matrix of the collision operator Ĉ.Secondly, the presented operator method could be adapted to other expansions of the single particle phase-space density function as well.(B2) Having derived an expression for the real matrix representation for an arbitrary operator, we can apply it to the angular momentum operators, the direction operators, the collision operator and the rotation operators.We note that O  contains complex terms and, thus, we have to check if the resulting representation matrices corresponding to the above list of operators are real.