{'fraction_non_alphanumeric': 0.03290570305988808, 'fraction_numerical': 0.01077509227288963, 'mean_word_length': 4.39226386325335, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Generalized Collisional Fluid Theory for Multi-Component, Multi-Temperature Plasma Using The Linearized Boltzmann Collision Operator for Scrape-Off Layer/Edge Applications

Grad's method is used on the linearized Boltzmann collision operator to derive the most general expressions for the collision coefficients for a multi-component, multi-temperature plasma up to rank-2. In doing so, the collision coefficients then get expressed as series sum of pure coefficients of temperature and mass ratios multiplied by the cross-section dependent Chapman-Cowling integrals. These collisional coefficients are compared to previously obtained coefficients by Zhdanov et al [Zhdanov V.M., Transport processes in multi-component plasma, Taylor and Francis (2002)] for 13N-moment multi-temperature scheme. First, the differences in coefficients are compared directly, and then the differences in first approximation to viscosity and friction force are compared. For the 13N-moment multi-temperature coefficients, it is found that they behave reasonably similarly for small temperature differences, but display substantial differences in the coefficients when the temperature differences are high, both for the coefficients and for viscosity and friction force values. Furthermore, the obtained coefficients are compared to the 21N-moment single-temperature approximation provided by Zhdanov et al, and it is seen that the differences are higher than the 13N-moment multi-temperature coefficients, and have substantial differences even in the vicinity of equal temperatures, especially for the viscosity and friction force calculations.

Power exhaust is a key challenge in next step fusion devices. Reducing the peak heat fluxes on the plasma facing components to tolerable levels, to a large extent, relies on impurity radiation in the boundary layer of the tokamak. The impurity radiation pattern in turn depends on plasma transport, both parallel and perpendicular to the magnetic field. This problem is generally addressed by solving plasma fluid models coupled to kinetic neutrals (e.g. as with the fluid code Soledge2d-EIRENE 1 ). Generally, the collision terms for the fluid dynamical equations are obtained by averaging over the kinetic equation with the constants of motion.
The moment-averaged collisional term can be determined by the use of different forms of ansatz for the distribution functions entering the collision operator of the kinetic equation. The two major ansatz used are namely, the Chapman-Enskog ansatz 2 and Grad's Hermite polynomial ansatz 3 . The Chapman-Enskog method involves decomposing the distribution function in terms of a small parameter expansion, and forming a momentaveraged hierarchy of equations at each order each with own collisional contribution. This method has been the most dominant so far, owing to its quick convergence. Grad's Hermite polynomial ansatz, on the other hand, involves decomposing the distribution function in terms of a series orthogonal tensorial polynomials, leading to a hierarchy of fluid equations for each order of the Hermite polynomial. These polynomials demonstrate mathematical properties which are somewhat easier to manipulate algebraically, however they have no clear rule for convergence. Often the convergence is either checked through brute force methods 4 , or by direct comparison of terms relevant to the physics in consideration 5,6 . Depending on the complexity of the collision operator, the treatment of the collision terms may get quite cumbersome.
Recently, in plasma physics oriented towards nuclear fusion, there has been a resurgence in the use of Grad's Hermite polynomial ansatz for use in calculating the collisional term, probably because of the improvement of algebraic techniques and availability of computer algebra systems. In recent plasma physics, it has been used to calculate the moments of the kinetic equation and the Landau collision operator 7 expressed in terms of Rosenbluth 8 potentials both for the linearized collision operator 9 and fully non-linear collision operator 10 , and its extension to magnetized plasmas 11 . Extension of moments of the Landau collision operator to strong flow cases has also been done 12 , and some closed-form analytic expressions for the involved integrals in the operator were also formulated 13 . The Hermite polynomial ansatz expressed in terms of the product of Laguerre polynomials and irreducible monomial, has also been used to formulate and study drift-kinetic models 14 and gyrokinetic models 15 for the scrape-off layer 16 , and also has been used to formulate a linear theory of electron plasma waves 17 .
However, the Landau collision operator is only valid for warm plasmas where the weak coupling conditions apply 18 . The generalization of these operators to different plasma regimes, for example in trying to account for shielding effects as in the Balescu-Lenard collision operator, leads to an increase in mathematical sophistication and a corresponding difficulty in solving 19 . The alternative is to use the much simpler Boltzmann collision operator. Given that the Landau operator can be thought of the Boltzmann operator under weak coupling limits 18 , we can expect that the use of the Boltzmann collision operator with the shielded Coulomb potential should provide quantitatively similar effects to the Landau or Balescu-Lenard collision operator as long as the plasma is in local thermodynamic equilibrium 20 and does not exhibit large scale fluctuations 21 . The added advantage of the Boltzmann collision operator is its use of explicit collision cross-sections. Any coefficients de-rived in this manner would have the advantage of being applicable for a wide variety of gas and plasma dynamics merely by using the relevant cross-section for the system in question 22 , for example for ion-ion [23][24][25] , ionneutral 26 , neutral-neutral 2,27-31 , and charge exchange collisions 32,33 .
In the scope of this article, we focus on the derivation of fluid collision coefficients from the linearized Boltzmann collision operator. In the past, two sets of collisional terms, one for any temperature range and the other for temperature range close to common plasma temperature, have been derived and provided 6,34,35 . However, an explicit derivation process was not provided for the values of the collision coefficients. Therefore, in order to verify the accuracy of the coefficients provided, we firstly rederive the collision operator in terms of partial bracket integrals, and derive the exact values of the partial bracket integrals 2,36 . We provide, for the first time, expressions for calculating the general collisional terms up to rank-2, in a manner that can be implemented efficiently in modern computer algebra systems. We also explicitly provide the range of validity of our and the aforementioned coefficients, and clearly delineate the underlying assumptions. This would be useful in clearly defining the simulation parameter range for a number of code packages which have implemented certain versions of the previous collisional coefficients. For example, the previous two sets of coefficients, taken from Ref. 6 , have been implemented in B2/SOLPS 37-39 , EDGE2D 40 , and more recently in Soledge2d-EIRENE 41 , which solves an energy equation for each species (while other codes mentioned solve only one total energy equation).
The article is organized as follows. We first provide a small introduction to the moment-averaged Boltzmann kinetic equation and the corresponding Boltzmann collision operator in Sec.I. We demonstrate its conservation properties in the process. Then, we introduce the Hermite polynomial ansatz and Grad's method, including the expression of the ansatz as a product of the Sonine polynomials and the irreducible tensorial monomial in Sec.II. In Sec.III, we present the derivation of the most general collision operator in terms of the partial bracket integrals. In Sec.IV, we provide the general forms of the partial bracket integrals (with the full derivations in Appendix A), alongwith the formulation of the cross-section integrals. Then we compare our obtained expressions for the collision coefficients with previously derived expressions in Sec.V for collisions of various species relevant to fusion, ranging from a deuterium-tritium plasma collisions to heavy impurity collisions such as with tungsten. Finally, in Sec.VI, we compare calculations of approximate values of physically intuitive quantities such as viscosity and friction force, and provide recommendations on the range of validity for the sets of coefficients.

I. THE BOLTZMANN EQUATION FOR MULTI-SPECIES PLASMA
The Boltzmann equation for the distribution function for species α, f α in the frame of the peculiar velocity of species α, c α = v α − u, is given by (1) where the common plasma flow velocity u is given by where ρ represents the mass density.The d/dt represents full time derivative given by d/dt = ∂/∂t + u.∇, and where the force term F α and du/dt are combined to write the relative force in the moving frame given by F * α = F α − m α du/dt.
The LHS is referred to as the "free-streaming term", and the RHS is the collisional contribution between species α and every other species of the system. The general "gain-loss" type Boltzmann collisional RHS is given by, where α, 1β refer to species of the two particles colliding, subscript ′ refers to properties after the collision, g is the relative velocity between the colliding particles, σ αβ is the collision cross section, and Ω is the solid angle in which the collision occurs. For the specific case of multispecies system, it takes the form where the distribution functions for each species are only dependent on the velocity of the species. Such a form is valid for elastic collisions. Now, for any quantity ψ α depending purely on species α, one can average over Eq. 1 which attains the following form d dt n α ψ α + n α ψ α ∇.u + ∇.(n α ψ α c α ) where n α is the number density of the species α, and where, For elastic collisions, the moment averaged collision operator can be transformed into since the distribution functions for any given species are purely a function of the species peculiar velocity. By its form, the averaged Boltzmann operator is meant to conserve mass, and by the choice of velocities, it is meant to conserve energy and momentum. One can notice in the averaged collision operator of the form Eq. (5), choosing ψ α = m α leads to a strict conservation of mass for R αβ for the averaged kinetic equation of each species. Hence, mass is strictly conserved. However, momentum and energy conservation can only be demonstrated over the sum of the averaged right hand sides of the kinetic equations for all species, i.e. α,β R αβ = 0. In order to demonstrate this, it would be sufficient to show that R αβ + R βα = 0. This is as follows One can notice that for momentum and energy, the term ψ ′ α +ψ ′ β −ψ α −ψ β vanishes as a result of the elastic nature of the collisions. Therefore, we now see that the Boltzmann collision operator is constructed to conserve mass, energy and momentum (and any linear combination of the three, for that matter), for any arbitrary form of the distribution function f . This fact is also useful to check the validity of the solutions obtained in the succeeding sections, acting as another check against calculation errors. We now proceed to choosing an ansatz in order to expand the collisional term.

II. SONINE-HERMITE POLYNOMIAL ANSATZ AND GRAD'S METHOD
In this section, we describe the modification of Grad's method 3 as used by Zhdanov 6 in his previous papers. In the ansatz for the solution of the Boltzmann equation, it is assumed that the solution f α is already near thermodynamic equilibrium for species α, f (0) α as follows where γ α = m α /kT α . In order to solve the Boltzmann equation (1), Zhdanov and Yushmanov choose an ansatz of the form Here, S n m+1/2 are the Sonine polynomials, given by, where the first few S n m+1/2 are Further, P (m) are the irreducible projection of the tensorial monomial c m α = c αr1 . . . c αrm , derived by the following recurrence relation with P (0) = 1. This expression is a sum of an outer product and the gradient with respect to the first-rank tensorial monomial γ 1/2 α c α . The first few P m+1 (γ 1/2 α c α ) are given by As one can observe, each P (m) is a rank-m irreducible tensor. We don't need to calculate any more than rank-2 for the scope of the current work. The constant τ mn arises as a result of internal contractions between b mn α and G mn α , and is given by The forms mentioned in Refs. 3 and 6 are cosmetically different because of the choice to use full factorial representations of functions and because of summing over full indices rather than over half-indices, but they are exactly the same. The coefficients b mn α are calculated as Some values of G mn α are as follows and the corresponding b mn α are given by Here, b 00 α , b 01 α and b 10 α represent the intuitive hydrodynamical moments density, diffusion velocity and temperature ρ α , w α = u α − u, and T α . the higher moments b 11 α and b 20 α represent the thermodynamically privileged moments (as per Balescu's nomenclature 42 ), the heat flux h α and the divergence free pressure-stress tensor π α , which are privileged because they contribute to the entropy. The higher-order moments b 12 α and b 21 α are nonprivileged moments r α and σ α , which do not have a clear physical meaning, however which may contribute to the accuracy of moment equations in terms of representing the Boltzmann equation. As one can notice, these are all moments of ranks less than 2. Moments of rank-0 are scalar, like density ρ α and temperature T α , are constitute N variables each. Moments rank-1 are vectorial moments, like momentum m α w α and heat-flux h α , and contribute 3N variables each. Moments of rank-2 are tensorial in nature, like the stress-tensor π α and σ α , and contribute 5N variables each (and not 9N , since they are symmetric and traceless). In principle, one can construct a 5N -moment system with just the hydrodynamical moments, a 13N -system with including the thermodynamically privileged moments, and a 21N -system including r α and σ α .
Furthermore, the Sonine-Hermite polynomials are chosen by Eq. (8) so as to form the moments in the most physically intuitive manner possible. They are related to the regular Sonine-Hermite polynomials in Ref. 3

, as follows
which then reduce to the Hermite polynomials defined in Refs. 43 and 44 by the following relation One limitation to note, is that the choice of a zero function f (0) α , defined at the common flow of the plasma, is valid for any difference of temperatures among different species, but limits the solution to requiring flow velocities of all species being approximately the same when the number of moments retained is finite. In principle, if one retained infinite moments, then the solution space provided by Eq. (7) would be the same as one provided by choosing an f (0) α defined at the individual species' flow velocity 45 . However, since we truncate this series at a very low number of moments, the solution spaces no longer overlap. However, the assumption of flows being close to each other is valid for SOL/edge of tokamaks, since the exit velocities of all species are close to the sound speed c s 46 . In general, one must keep in mind the general ordering of the diffusion velocities as follows One can now introduce ψ α = G mn α in the averaged Boltzmann equation (3), in order to compute the moments, obtaining an infinite hierarchy of transport equations.
Generally, in the LHS of the hierarchy of balance equations for the moments obtained in this manner, the k th moment equation contains both the (k−1) th and (k+1) th moments. Hence, one has only k equations for k + 1 variables. In order to provide a closure, in Grad's method, one truncates at a moment k, such that moments higher than k are calculated using the expansion for f α truncated at the k th term, which approximates the higher moments in terms of the lower ones. This closes the set of equations obtained. Illustration of Grad's closure and also Zhdanov closure are out of the scope of the current article, and will be demonstrated in an upcoming article. In this article, we illustrate the development and solution of the RHS of the equation, i.e. the moment averaged Boltzmann collisional operator, and compare it to previously obtained values by Zhdanov et al 6,34,35 .

III. DERIVATION OF THE RIGHT HAND SIDE OF THE BOLTZMANN EQUATION
In the Boltzmann collision integral J αβ , it is possible to choose a distribution function which takes the form , which essentially represents the ansatz as a perturbed Maxwellian. The moment-averaged collision operator Eq. (5) can be written as on neglecting the ΦΦ terms. This linearization of the collision operator makes it bilinear, i.e. they satisfy the relations and correspondingly the moment-average R(ψ, f α , f β ) is trilinear.
This allows us to decompose the moment-average into sums of smaller terms, which is useful analytically. This also is similar to the properties exhibited by some other linearized operators such as the linearized Landau operator 42,47 . On substituting the Sonine-Hermite polynomial ansatz from Eq. (7) for the distribution functions f , and set the moment ψ = G mn α from Eq. (8), we obtain where R mnkl αβ represents the part of R mn αβ averaging over the kl term of the ansatz Eq. (7). Noting that G kl b kl is an inner product, we now substitute the definition of G mn , and use the following integral identity 9,48 where W is symmetric and traceless tensor of rank k not a function of v. Furthermore, we define a "bracket" integrals of the following form through which we can contract over index k in R mnkl αβ and write it as R mnl αβ , such that where A mnl αβ , B mnl αβ and C mnl αβ are given by where Q mnl αβ = (−1) n+l 2 2l+m (2m)!(m + l)!n! (m!) 2 (2m + 2l + 1)! n α n β . (16) This moment-averaged collision operator is valid for any difference of masses or temperatures of the colliding species.
Having now derived our expressions for the momentaveraged collision operator, we now illustrate a slightly different expression derived originally in the appendix of Ref. 6 by Zhdanov, which brings with it some additional assumptions. The collision operator is linearized by choosing a distribution function defined at the common temperature of the plasma T = α ρ α T α / α ρ α , such that f on neglecting the squared ΦΦ terms. However, in such a form of the collision operator, Zhdanov et al assume that the temperatures of each species is close to the common temperature, i.e. |T − T α | ≪ T . One could term this as the linearized Boltzmann collision integral under quasi thermodynamic equilibrium conditions. Such a scheme is well-suited when the masses of the colliding species are of the same order, thus both species possessing a collisional relaxation timescale of the same order as well, leading to their distribution functions being close enough to a Maxwellian defined at the common temperature T . However, one needs to be careful when the masses of the species are at different orders, leading to different relaxation timescales, such as in the case of heavy impurities colliding with the plasma fuel species, where relaxation timescales may be different. On integrating over a moment ψ α = G mn α to obtain the moment averaged collision integral using the ansatz Eqs. (7) for the distribution function, and on noting that b ′kl α = b kl α for a given species α by construction since the form of f α remains the same, we can obtain where A mnl αβ and B mnl αβ are given by where Q mnl αβ has the same expression as Eq. (16). This is also the collision operator employed for deriving the 21N -moment single-temperature collision coefficients in Ref. 35. Note, however, that the A mnl αβ , B mnl αβ here correspond to the A mln αβ , B mln αβ of Eq. (15). The limit of the moment averaged collision operator (17) is that it is only valid for quasi thermodynamic equilibrium conditions with the species temperatures being close to the plasma common temperature. However, the one in Eq. (14) has no such assumption. Therefore, each individual moment averaged term of Eq. (17) will be less accurate than Eq. (14) for increasing difference in the temperatures of the colliding species.
It is worth discussing what the assumption of b ′kl = b kl means. Since b kl enter the distribution functions, this assumption essentially implies no change in the physical quantity in the pre-collision and post-collision distribution functions. For example b 01 α = ρ α w α is the diffusion velocity of the species α, and assuming that b 01 α = b ′01 α implies that the diffusion velocities of the pre-collision and post-collision distributions remain the same for the same species α, not changing in the timescale of the collision. This would also ensure that the distribution function f has the same form for a given species α for all four pre-collision and post-collision distributions, as necessitated by the Boltzmann equation. This is a reasonable assumption when the duration of the collision is very small as in the case of short-range forces, e.g. in case of rigid-sphere collisions where the collision only lasts an instant, or for weakly-coupled long-range forces such as Coulomb potential, where again the small-angle collision duration is very small. However, for long-range interaction potentials with strong coupling, one has to be careful that the collision duration is much smaller than the timescale of the system evolution, and in general for such an interaction potential the Boltzmann collision operator is anyhow not appropriate.
To use either expression of the moment-averaged collision operator, we need to calculate these bracket integrals for the required (m, n) values in order to derive the desired forces. Because of Eq. (12), the only brackets that survive in the moment-averaged linearized collision operator are the ones possessing the same rank, which significantly reduces the number of terms one must calculate. Certain methods for deriving these bracket integrals are provided in Refs. 2 , 5 and 36 for m = 1, 2, as they are well suited enough within the scope of the 21Nmoment scheme. In the next section we indicate the general expressions for the bracket integrals obtained by following these methods for a case of species α, β possessing different masses and different temperatures.

IV. GENERAL EXPRESSIONS FOR THE BRACKET INTEGRALS
The general expressions for the rank-m bracket integrals take the following general forms The coefficients A pqrl,m αβ , A pqrl,m αα are functions of mass and temperature ratios of the species α and β. The terms Ω lr αβ are the effective cross-section moment integrals of Chapman and Cowling (henceforth referred to as the "Chapman-Cowling integrals") which are dependent on the potential of interaction between species α and β. The factor k αβ is an arbitrary function of the masses and temperatures of the colliding species, which can be chosen freely as long as k αβ > 0. The choice also affects the forms of A pqrl,m αβ ,A pqrl,m αα and Ω lr αβ which depend on k αβ individually, such that, in principle, the overall bracket integral values do not depend on k αβ . The exact derivations the generalized coefficients A pqrl,m αβ and A pqrl,m αα for different bracket integrals up to rank-2, with all steps supplied for verification purposes, are provided in Appendix A. For the purpose of our work, up to rank-2 suffices, as in the context of the Boltzmann collision operator, only the moments up to rank-2 have been considered in previous works. The solution for higherorder bracket integrals is out of the scope of the current article, and will be reserved for a future manuscript.

The
"reversed" bracket integrals S p 3/2 (W 2 β )P (m) (W β ), S q 3/2 (W 2 α )P (m) (W α ) required for the second form of the collision operator, and S p 3/2 (W 2 β )P (m) (W β ), S q 3/2 (W 2 β )P (m) (W β ), required for verifying conservation properties, can be obtained by transforming α ⇋ β in the above expressions provided. Since the two types bracket integrals are calculated independently of each other, any mistake in calculating them would lead to non-conservation of momentum or energy. Hence, demonstrating conservation of momentum and energy with the obtained quantities is an adequate testament to their accuracy (see Appendix B). Furthermore, the expressions for the bracket integrals are very amenable to being implemented in computer algebra systems, as they are just composed of sums and products which any computer algebra system should be able to perform. We implement the expressions for the bracket integrals and the moment-averaged collision term in Mathematica 49 .
The Chapman-Cowling integral is written in the following form for our case where ζ = d 1/2 αβ g and where the factor d αβ is related to k αβ by the following relation Choosing k αβ fixes d αβ and vice versa, and note therefore that Ω lr αβ is a functional of d αβ and the effective cross-section φ (l) αβ . Essentially we are now left with calculating the effective cross sections φ (l) αβ , which in turn depend on the physics of the particle-potential interaction. Therefore, the choice of potential is crucial to calculating these effective cross sections accurately. We essentially have a choice between pure Coulomb potential or the shielded Coulomb potential, as these are the only ones that apply to fully ionized plasmas in the fusion domain. However, there are always integrability and convergence issues with the potentials used in these calculations. For example, the pure Coulomb potential diverges at the limit of low collision angles with high impact parameters which constitute the majority of collisions in a hot plasma. This is usually mitigated by choosing a cutoff for the integral at the Debye length radius, which then leads to a converged integral. For the shielded Coulomb potential, for high energy and low angle interactions, it may lead to an issue where some forms of the integrals obtained do not converge at high impact parameter limits. Some physical approximations, such as ignoring large angle collisions for the shielded Coulomb potential part of the integration, are used to express the integral approximately in forms that converge. From a modelling point of view, it is worth keeping in mind that these two potentials will offer slightly different collisional coefficients, which can provide a range of values which could be useful for comparison with experiments.
In general, the integration of the effective crosssections with the shielded Coulomb potential is not so simple, and often has to be manually done for different values of l, with unique approximations applied at each value. However, it is possible to find the asymptotic values of the cross-sections through a perturbation method 23 .
For our case, the asymptotic form of Ωintegral for the shielded Coulomb potential looks like where where T is the plasma common temperature given by where ∆ αβ is the mean distance of closest approach (also called the "particle diameter"), λ D is the Debye length given by and γ is the Euler-Mascheroni constant. This formula has the advantage of being easily implementable in computer algebra softwares. However, it is worth noting that this expression leads to the Chapman-Cowling integral for the shielded Coulomb potential (and hence also the bracket integrals) evaluating to different values depending on the chosen value of d αβ . This arises from two assumptions in the calculation of the cross-section from the shielded Coulomb potential in Refs. 23 and 24. The first being the non-dimensionalizing of the impact parameter with respect to the particle diameter and the Debye length, which then vanishes under the asymptotic lower limit of the collision integral, leaving the final result the same for any choice of the Debye length or the particle diameter. The second is that the model of collisions is essentially a collision model of one particle being deflected by one potential. The cross section integrals φ (l) , found in this manner depend only on the relative velocity g other than the potential of interaction. However, the relative velocity g in the Chapman-Cowling is scaled by a factor d 1/2 αβ , and hence the value of φ (l) will also need to be transformed to the scaled value of ζ = d 1/2 αβ γ, leading to the additional factors of 2kT d αβ /µ αβ seen in Eq. (21). For our coefficients, we choose d αβ in such a manner that d αβ = µ αβ /2kT , in order to agree with the two aforementioned calculations of the cross-sections for the shielded Coulomb potential in Refs. 23 and 24. Having d αβ in such a manner also allows us to immediately use cross-section values from numerical and empirical databases such as AMJUEL 50 , LXCat 51 and ADAS 52 , because the cross-sections in these databases are often expressed as a polynomial of µ αβ g 2 /2kT . Some attention also needs to be paid to the chosen value of ∆ αβ . Following Liboff 24 , Kihara et al 23 where T is the common temperature of the plasma. Such a value of ∆ αβ is also termed the "particle diameter" in the literature (where π∆ 2 αβ would be the collision crosssection for the rigid-sphere collision case 2 ). Following Zhdanov 6 , one could choose it as the average inverse impact parameter 1/b 0 −1 over the distributions as follows where γ αβ is given by γ αβ = γ α γ β /(γ α + γ β ). (which for equal temperatures makes the logarithmic term in the expression equal to the Coulomb logarithm ln Λ αβ ). Whichever is chosen has to be used consistently. In Ref. 6 by Zhdanov, the cross-section has been treated by using Eq. (23) for the Coulomb logarithm, but Eq. (22) for the mean distance of closest approach outside the logarithm. This is essentially the same approximation as in Rosenbluth et al 8 , but which has non-negligible effect on the Ω lr αβ for high values of (l, r) (since neither A n nor C n converge for large n, and can become larger than the log term for a large enough value of n). This still remains a reasonable approximation however, because we are only concerned with a small number of lower-order moments.
One can also, following the interpolated formula provided above, extrapolate results obtained by following Liboff's procedure 24 for results provided by Bonnefoi 53 , and one will find the same expression as Eq. (21), but with −2γ in place of the −2 ln γ. The cause of this is that, in Liboff's procedure, the general angle of deflection of the collision takes the form of the first-order modified Bessel function of the second kind K 1 , and one must use approximations on the order of the energy of the interaction to express the integral at any order as the square of K 1 , as every other power of K 1 is nonintegrable in the given limits. This leads to the integral having a correction of A r − γ instead of that in Eq. (21).
One feature to note of these shielded Coulomb potential cross-sections is that they ignore any difference in shielding for the attractive and repulsive potential cases, because the particle diameter is assumed to be positive for both attractive and repulsive cases. Depending on the kind of plasma, this may be relevant 22,54 . However, for the case of fusion plasmas, this may safely be neglected, given that the logarithmic term is dominant for hot plasmas.
We also mention for the sake of completeness that the original formula for the Chapman-Cowling integrals used by Zhdanov for the values of Ω lr αβ in Ref. 6 found using the Coulomb potential with the Debye length cutoff is given by where the common plasma temperature T refers to T = (1/ρ) α ρ α T α . This follows from Zhdanov's forms of the Ω-integrals, when the temperatures of the components are close to each other and when the Debye length λ D is much smaller than the inverse of the average inversed impact parameter, i.e.
and where Λ αβ is the Coulomb logarithm given by This follows the general result in Eq. (21) to the order of the logarithmic term. We will calculate our cross-sections with Eq. (21) with Eq. (22) and with d αβ = µ αβ /(2kT ). With expressions for the coefficients A pqrl αβ/αα and the cross-sections Ω lr αβ,sh , we can now proceed to calculate the collisional coefficients and compare them to the ones found in existing literature.

V. RANGE OF VALIDITY OF ZHDANOV'S VALUES
Zhdanov et al have previously derived two sets of coefficients for the values of the collisional coefficients. The first set was derived for a multi-temperature, multicomponent plasma without any explicit assumptions on the temperatures of the species (see Appendix C). These 13N -moments multi-temperature coefficients were also calculated using the linearized Boltzmann operator. The derivation method for a few lower order moments is given in Refs. 6 and 55. However, to the best of our knowledge, no explicit general derivation scheme method was provided, making it difficult to verify some of the cumbersome higher-order moments, and generate higher-order ones, in case we need larger number of moments. These coefficients generally take the form where K m,pqrl αβ,α is a term depending solely on the masses, temperatures and densities of the species, Ω rl αβ is the Chapman-Cowling effective cross-section moment defined by Eq. (C28), and Θ αβ is given by We can verify that on choosing d αβ = γ αβ /2, our coefficients equal the ones provided by Zhdanov (see Appendix C), thus verifying the coefficients provided in known literature for the chosen value of d αβ . We can now also generate higher-order moments of this multitemperature form if required. The second set of collision coefficients, provided by Zhdanov in Refs. 6 and 35, are derived using the form of the collision operator given by Eq. (17), for 21N -moments (see Appendix D). They take the general form where L m,pqrl αβ,α is a term depending on the masses, temperatures and densities of the species. The bracket integrals are evaluated at the common plasma temperature T = α ρ α T α / α ρ α , however. There still remains, however, an individual species temperature dependence from the terms multiplying the bracket integral (See Exs. (18)). The values of Ω rl αβ , however, are proportional to the approximate formula Eq. (D1), which is accurate to the order of the Coulomb logarithm.
When comparing our expressions derived with d αβ = µ αβ /(2kT ) with the expressions Eqs. (25) and Eqs. (26), we can expect the 13N -moment multi-temperature coefficients to agree in the case of equal temperatures and the 21N -moment coefficients to modestly agree, and for both have a range of reasonable agreement in the vicinity of equal temperatures. To understand where the coefficients begin to diverge, we consider some physical situations relevant to SOL/edge physics, i.e. 1. a three component plasma intended to be the fusion fuel, with electrons, deuterium (D) and tritium (T), the three being at comparable densities, as D-T fusion is planned to be used in current and future high-Q campaigns, 2. a three component plasma with light impurities at significant fraction (10%) of the main fuel species, i.e. electrons, hydrogen (H) and carbon (C), with the carbon in the plasma originating from facing plasma components made of graphite, 3. a three component plasma with injected mid-weight impurities with densities at a small fraction (1%) of the fuel species density, electrons at a small fraction, hydrogen and argon (Ar), often used for controlled experimentation with impurities or for other purposes, and finally 4. a three component plasma with a heavy impurity at trace levels (0.001%), i.e. electrons, ions and tungsten (W), where the tungsten usually originates from the walls and divertors made of tungsten. Such a choice of scenarios will help us scan over operationally relevant mass ratios and density ratios, allowing us to focus on the effect of the temperature ratio. The values of masses, charges, densities and temperatures we choose can be found in Table (I). The temperatures are chosen so as to provide a range of temperature ratio spanning 0.1 − 2.
Firstly, we begin by looking at coefficients for the D-T case with a physical significance that can be intuitively understood, e.g. the friction force, governed by A 100 ij , B 100 ij , the thermal gradient force, which is governed by A 110 ij , B 110 ij , the energy exchange term given by C 010 ij . We can immediately notice in Fig. (1), that generally in the vicinity of equal temperatures T D /T T = 1, the curves for all coefficients tend to follow each other quite closely. In particular, the curves for the 13N -moment multi-temperature coefficients follow our obtained values much closer than the 21N -moment single-temperature ones. However, they deviate quite significantly going away from equal temperatures. Based on these observations, we state that the 13Nmoment multi-temperature coefficients agree more with ours in the vicinity of equal temperature than the 21Nmoment single-temperature ones. However, to better recommend a range of validity and quantitatively characterize the deviations, we proceed to plot the percentage differences in the coefficients, defined as the absolute value percentage difference of Zhdanov's two sets of coefficients with respect to our coefficients Eqs. (14).
The percentage differences in the 13N -moment multitemperature coefficients are plotted in Fig. (2). We can notice exact agreement at equal temperatures, and that the differences in the coefficients are very low in the vicinity of equal temperatures. However, they seem to deviate rapidly as the temperature ratio decreases below 0.5. In particular, the coefficients related to the heat flux transmission A 110 αβ , B 110 αβ seem to deviate very quickly. This would indicate that at significant temperature differences, the representation of heat flux gains more importance. In Table (II), we showcase the maximum differences in the temperature ratio range of 0.8 − 1.2. We can see that all differences for all four physical cases are less than 22%, with most coefficients having differences less than 6%. Furthermore, we can notice that the heavier the impurity becomes, the lower the differences are. This indicates that the 13N -moment multi-temperature coefficients may be more suitable for simulation the heavier the impurity species being simulated.
We show the differences in the single-temperature 21N -moment single-temperature coefficients for the D-T case in Figs. (3), ((4)) and (5). One can notice in these figures that the differences in the coefficients are significantly higher than those of the 13N -moment multitemperature case. Furthermore, the percentage differences in higher order moments, e.g. A/B 11n αβ , A/B 12n αβ and A/B 21n αβ are significantly higher than those of the lower order moment ones. The same trends are observed in the cases of carbon, argon and tungsten as well. In order to compare the differences, in Table (III), we show the maximum difference in the temperature ratio range of 0.8 − 1.2. We notice the same trend as earlier for the lower order coefficients, in the decrease in the percentage In the legends for the plots, first A mnl αβ /B mnl αβ refer to our coefficients, the "13N " indicates Zhdanov's multi-temperature collisional coefficients and similarly the "21N " refers to Zhdanov's single-temperature collisional coefficients.  Plots of percentage differences in multi-temperature coefficients calculated by Eqs. (25) and ours calculated from Eq. 14. The top plot showcases the difference for the rank-0 coefficient, the two middle plots shows the differences for the rank-1 coefficients, and the bottom plot shows the differences for the rank-2 coefficients. The spike near unity temperature ratio is because of our collision coefficient changing signs, but the error being the absolute value of the relative difference. differences the heavier the impurity species gets. However, we also notice that the difference for the higherorder A/B 11n αβ , A/B 12n αβ and A/B 21n αβ moments can be up to three orders of magnitude higher than the lower order ones. Thus, the 21N -moment single-temperature coefficients cannot be recommended for simulation purposes with significant temperature differences, compared to the 13N -moment multi-temperature ones. The top and the middle plots showcase the differences for the rank-1 coefficients, and the bottom plot shows the differences for the rank-2 coefficients.

VI. EFFECT OF COEFFICIENTS ON VISCOSITY AND FRICTION FORCE CALCULATIONS
We have up until now clearly noted the difference in coefficients numerically. However, it may also be instructive to study some intuitive physical quantities such as viscosity and heat flux. The process of obtaining values of viscosity and heat-flux would essentially close the 13N -moment system of equations given by (3), by eliminating the higher-order moments, h and π in this case, in terms of lower order ones.
In order to obtain the closed set of equations, we neglect any electric and magnetic fields, and follow the procedure in Ref. 56, restricting the collisional terms The top and the middle plots showcase the differences for the rank-1 coefficients, and the bottom plot shows the differences for the rank-2 coefficients.
to the 13N -moment multi-temperature approximation. Clearly, this approximation will yield viscosity and heatflux purely inertial in origin. In addition to neglecting the fields, we also assume that the higher order moments evolve much slower in time and have much smaller gradients than the lower order plasma dynamical moments. In terms of the Knudsen number Kn ∼ λ/L ∼ τ /T where λ refers to the mean free path and L the scale length of the system in question, and equivalently τ is the collision frequency and T the timescale of the evolution of the system, we see, therefore that the plasma dynamical moments are of the order of the Kn, and ignore every term that is of a higher order than Kn. This implies that we neglect all time and space gradients of higher order moments in the LHS of the moment-averaged Boltzmann equation (3). Hence, the only terms that survive are the ones proportional to dG mn α /dt (through the T α dependence of γ α ) for rank-1 quantities, and the ones proportional to ∂u r /∂x s for the rank-2 quantities (as the c αs derivative reduces the order of the moment).
Thus, the reduced evolution equation, now the steadystate equation, for the stress tensor π α for the species α becomes where ε is given by This is the usual form by which one represents viscous forces as the viscosity η multiplied to a traceless strain rate tensor ε 57 . We can re-write this equation in the following form where the sum γ is over all species including α, and where For N species in the plasma, the set of N equations corresponding to Eq. (27) for all species can then be compactly written in a matricial form where the B * 200 is an N × N matrix whose elements are given by where ω is an arbitrarily chosen n th species, and where Π and P is the column matrix with all the values of π γ and 2p γ respectively. On inverting the equation, and comparing it to the classical form of stress tensor π γ = −2η γ ε, we can find the column matrix of the viscosity E given by where the elements of E are the partial viscosity of each species η γ . One can find the total viscosity η = γ η γ by the formula η = T r(E T U ).
Similarly, the reduced heat-flux h α equation is given by which can then be written in a matrix form as One can notice that approximating the heat-fluxes in this manner, reduces them to a linear combination of temperature gradients and flows. On substituting this value of the heat flux, the RHS of the momentum balance equation will depend solely on only two terms, one proportional to the flows and the other proportional to the temperature gradients, hence recovering the familiar form of the collision term, one with friction force dependent on the flow difference and the other a thermal gradient force, dependent on the difference of temperature gradients 46,58 . The heat-flux term then adds to the existing friction force term and augments it in the following manner and the thermal gradient force remains where [. . . ] γ indicates the element of the column matrix corresponding to the species γ. To estimate the augmentation of the friction force, we need to compare the additional term to the original coefficients A 100 αβ m α and B 100 αβ m β . One can notice, however, that the friction force term R 1 αβ,f ric now becomes quite complex, with friction among any two pairs of species depending on the flow velocity of all species. This makes any straightforward comparisons of the coefficients of the friction force cumbersome and prone to overinformation. Therefore, we only choose the coefficient of w α and w β from the heatflux contribution h α and h β and compare them to the  6. Plot of total viscosities for the Deuterium-Tritium plasma calculated using our coefficients, and Zhdanov's 13N -moment multi-temperature and 21N -moment singletemperature coefficients, plotted against the temperature ratio TD/TT . Notice the virtual overlap between the viscosities calculated from Zhdanov's 13N -moment coefficients and ours in the vicinity of equal temperature. In the legends for the plots, first η refers to our coefficients, the subscript "13N " indicates the viscosity calculated Zhdanov's multi-temperature collisional coefficients and similarly the subscript "21N " for Zhdanov's single-temperature collisional coefficients.
original coefficients corresponding to the original friction force. Furthermore, one can repeat the same procedure with a higher-order number of moments as has been done in Ref. 35, often referred to as the "Zhdanov closure" when applied to close a 21N -moment set of equations with the moments of the collision operator given by Eq. (17). It must be mentioned that the higher approximations to the viscosity, which would depend on the electric field and the magnetic field, can be decomposed into a form which is a linear combination of different viscosity contributions 34,58,59 . Comparison of such approximations with more nuanced field effects and higher number of moments is out of the scope of the current article, since a larger number of moments would imply needing the use of block matrices instead of regular matrices used here, and will be a part of our planned future work.
For the particular case of a Deuterium-Tritium plasma, one can notice the total viscosity η in Fig. (6) in cases of all three coefficients follow each other quite closely in the vicinity of equal temperatures, with the 13N -moment multi-temperature ones practically overlapping with our exact values. One can see in Tables (IV) and (V) that the differences in viscosity for the 13Nmoment multi-temperature case are significantly lower than those of the 21N -moment single-temperature case, and also that differences in viscosity generally seem to decrease with increasing mass ratio.
Furthermore, we can compare the viscosity values obtained with the prescription for parallel viscosity given by Braginskii (Ref. 58, page 229), where Braginskii provides the magnitude of viscosity by the expression η ∼ nkT τ αβ , where τ αβ is the mean time between collisions between species α and β, and where n = i n i . Values of τ αβ are provided by Braginskii for ions and electrons, but not for impurities. However, since Braginskii follows the Chapman- where the first approximation to the diffusion coefficient [D αβ ] 1 is given by On substituting the physical values from Table I, and the value of Ω 11 αβ from Eq. (21), one can find that the viscosity prescribed by Braginskii, η ∼ nkT τ αβ , falls within the same range of values as found in Fig. (6).
A similar trend can be partly observed in the augmentation of the friction force coefficients. In Fig. (7), the increase in friction force is indicated as a percentage of the original value. In general, with increasing mass ratio and decreasing densities, the augmentation in the term proportional to w α decreases and increases in the term proportional to w β . This is expected, as the dominant part of the heat-flux contribution to the friction force is bound to arise from the species with the larger density. In case of the 21N -moment single-temperature cases for the friction force augmentation for the term proportional to the background flow contribution, the differences seem to rise again for high mass ratios, as evidenced by the W-H case. The percentage differences on these computed physical quantities seem to have a minimum in between the Argon and Tungsten cases, which indicates that even the relatively smaller differences in the 21N -moment single-temperature case may contribute to significant difference in physical quantities of interest. On the basis of this, we recommend caution in using the 21N -moment single-temperature coefficients even for heavy impurities when temperature differences may be significant.

D-T C-H Ar-H W-H
Thus, from these numerical results we can conclude 1. that the coefficients calculated by our coefficients, Zhdanov's 13N -moment multi-temperature and 21Nmoment single-temperature expressions tend to follow each other quite closely in the vicinity of equal temperatures, 2. that the differences for the 21N -moment single-temperature coefficients are higher than those of the 13N -moment multi-temperature coefficients, 3. that the differences in coefficients become quite significant outside the vicinity of equal temperatures, however, that in general again, the differences in 13N -moment multitemperature coefficients are less than those of the 21Nmoment single-temperature case, 4. that the differences in higher-order moments become quite significant in case of the 21N -moment single-temperature coefficients, 5. that the differences in coefficients in general decrease with increasing mass ratio, and 6. that despite the agreement and trends being similar in physical quantities such as viscosity and augmentation in friction force, computed from the coefficients in the vicinity of equal temperatures, the values computed from the 21N -moment singletemperature coefficients may present significant differences as the mass ratio increases, especially in the augmentation of the friction force term.

VII. SUMMARY, CONCLUSIONS AND FUTURE WORK
In this article, we have generalized the calculation of moment averaged collisional coefficients for a multicomponent, multi-temperature plasma without making any assumptions on the masses or temperatures of the colliding species with the linearized Boltzmann collision operator for up to rank-2 tensorial moments. We started by taking an ansatz for the distribution function in terms of the Sonine polynomials and the irreducible tensorial Hermite monomials Eq. (7), and then expressing the moment in terms of these polynomials of different orders Eq. (8). On taking moments defined in this manner, and then averaging over the Boltzmann collision operator, we obtain a generalized moment-averaged collision term, which expresses itself in terms of partial bracket integrals in Eq. (14). We evaluate these partial bracket integrals analytically and derive general expressions up to rank-2. The collision operator found in this manner is valid for any range of masses and temperatures, but is restricted to the flow differences being much smaller than the order of the thermal velocities of the plasma. Furthermore, the collision term automatically preserves energy and momentum due to symmetry properties of the Boltzmann collision operator. These expressions for the collisional coefficients were found to be very amenable to being implemented in computer algebra systems, and in our case were implemented in Mathematica 49 , and were found to conserve mass, energy and momentum.
The collision coefficients are found to essentially be linear combinations of product of a term A pqrl αβ/αα depending purely on masses and temperatures, and another purely depending on the collisional cross-section Ω lr αβ , which in turn purely depends on the potential of interaction between the two colliding species and a factor d αβ . For the cross-sections, we chose a formula for the cross-section derived from asymptotic values of the form of the crosssection integral for the shielded Coulomb potential from Eq. (21). We calculated our set of coefficients choosing d αβ = µ αβ /(2kT ), as choosing this factor agrees with the calculations of effective cross-section integral for the shielded Coulomb potential from previous literature.
Previously, Zhdanov et al had derived two sets of expressions for the collision coefficients and cross-sections, the first being derived for a multi-temperature plasma with terms provided up to 13N -moments, and the second being derived at the common plasma temperature, with terms provided up to 21N -moments. The procedure of derivation for the former is not explicitly provided for all moments, and no expressions for the coefficients were provided for higher-order collisional moments.
However, with our calculation procedure for the bracket integrals, with the suitable factor d αβ = γ αβ /2, the provided moments were found to be accurately derived. The procedure for the latter 21N -moment single-temperature set of coefficients was explicitly provided, but comes at the cost of temperature differences being much smaller than the order of the plasma common temperature.
We then compared our expressions to the ones provided by Zhdanov for multiple cases of colliding species relevant to fusion. We find that in the vicinity of the temperatures being equal, all sets of coefficients agree very well, but they diverge away from equal temperatures. We further find that the differences in the 13Nmoment multi-temperature coefficients are much smaller than those of the 21N -moment single-temperature ones (See Appendices C and D). We also find that the differences in the coefficients decrease as the mass-ratio of the species increases. Furthermore, we use certain approximations in Knudsen number to obtain reduced equations for the stress-tensor and heat-flux, and use it to calculate the first inertial approximation to total viscosity and the augmentation of the friction force as contributed by the heat-flux. We find that while the differences mostly follow the same patterns as the coefficients, for high mass ratio and low density, as in the case of tungsten impurity in tokamak plasmas, the differences in these physical quantities become significant for the 21N -moment single-temperature case. On the basis of this, we caution against using the 21N -moment single-temperature coefficients where the temperature difference between species is significant. In the same vein, we find the 13N -moment multi-temperature coefficients are agree better than the 21N -moment single-temperature ones for small temperature differences, especially for modelling heavy-impurity transport, as for tungsten. For any significant temperature differences between species, we recommend using multi-temperature coefficients.
In the future, we plan to generalize the expressions for the linearized collision operator up to arbitrary rank-m tensorial moments. We also plan to study in detail linear parallel closure schemes associated with fluid equations derived from the moments under the low Knudsen number approximation. The ones with the aforementioned 21N -moment single-temperature scheme, derived by Zhdanov et al, often termed the Zhdanov closure, is one such scheme of closure particularly of interest to us, since it has already been implemented in many fluid codes of interest to the fusion community.

Appendix A: Derivation of the bracket integrals
We follow the method illustrated in Refs. 2 and 36 in order to derive these bracket integrals in their most general form. Certain transforms in Ref. 36 are taken from R.S. Devoto's thesis 60 . We proceed in a detailed manner to show all the calculations necessary to calculate the brackets.
There are several transforms one must make to render terms in the bracket integral tractable. The first transform is to move to the center-of-mass frame, where the center-of-mass velocity relative to the common flow U and relative velocity g are defined as where the superscript (′) represents the same corresponding expression for the primed variables, and consequently where µ αβ is the reduced mass given by µ αβ = m α m β /(m α + m β ). Now we define variable changes X andḡ where The Jacobian of the transform is given by These transforms will help transform the velocities into quantities that are easier to integrate. In order to deal with the Sonine polynomials, it is worth noting that the Sonine polynomial S n m (x) is defined to be coefficient of s n in the expansion of With this relation in mind, we now proceed to derive some bracket integrals. We describe the derivation of S p 3/2 (W 2 β )W β , S q 3/2 (W 2 α )W α in full detail, after which we provide condensed derivations for the other brackets with only the differing steps.
We have the bracket integral defined as follows In order to recover collision coefficients out of this form of a bracket integrals, the general strategy will be to use the property (A7) and absorb the Sonine polynomial terms into the exponential, creating a generating function Π αβ which contains the values of the bracket integrals for all (p, q). Then, the term left outside of the exponential would just be the term corresponding to the full contraction of the irreducible Hermite polynomials. Certain transformations will be performed to integrate this full contraction term over the exponential. However, since the Sonine polynomials were absorbed into the exponential, the remaining exponential after the integral will have to be expanded into a series to isolate the (p, q) coefficients. Once the integral is performed, the terms left can be expressed as linear combination of moments of the relative velocity g over the cross section σ αβ and a gaussian, Ω lr αβ . The generating function Π αβ in found to be in the following form with the help of Eq. (A7), where S = s/(1 − s) and T = t/(1 − t), and where the bracket integral of order (p, q) is the coefficient of s p t q in the expansion of Π αβ . Now we substitute in the expressions for the zero functions f and make the coordinate transform (c α , c β ) → (X,ḡ), and we get A14) which in the limit of a "no-collision event" becomes (ḡ, χ = 0) such that This helps us split the generating function first as an integral over X leaving H αβ purely a function of (ḡ, χ). Thus, Π αβ as Now defining where θ αβ = T α /T β . Then H αβ (ḡ, χ) can be transformed to We now make another variable transformation as follows, so as to write the exponential term purely in terms of the squares of quantitiesX andḡ, such that the Jacobian of this transform, (X,ḡ) → (X,ḡ), is 1, and further defining b αβ = a βα − M α1 M α2 a αβ (S * 2 + T * 2 + 2S * T * cos χ), (A21) we can write H αβ (ḡ, χ) as (A22) At this point, it now remains to derive and expression for W ′ β .W α . In order to do so, we introduce V α and V β and their primed counterparts such that W α and W ′ β can be expressed in terms ofX andḡ (′) These relations are be obtained by substituting Eqs. (A20) in Eqs. (A3). Therefore, the dot product Now, odd powers ofX will vanish in the integral for H αβ (ḡ, χ), so we only need to concern ourselves with expressing V ′ β .V α back in terms ofḡ. This calculates to which on substituting back into W ′ β .W α , and defining so as to simplify W ′ β .W α as Therefore, now H αβ (ḡ, χ), on substituting W ′ β .W α and integrating overX, becomes Now that H αβ (ḡ, χ) has been evaluated, it remains to express Π αβ in terms of s and t, so that the bracket integral can manifest itself as the coefficient of s p t q . For this purpose, a number of series expansions need to be performed. We begin by performing a Taylor expansion of exp(−b αβḡ 2 ) in the following form, by separating out an arbitrary factor of exp(−k αβḡ 2 ) pre-emptively so that quantities in the integral can be expressed later as moments of a gaussian, and therefore (A31 Now, on substituting the expression for b αβ , expanding a αβ , a βα , S * and T * from their definitions, and setting S = s/(1−s) and T = t/(1−t), H αβ (ḡ, χ) can be written as follows where the N 's and D's are known coefficients of the form Now, we need to express the numerator and denominator as series expansions in order explicitly obtain and group together the powers of s and t. We proceed to expand the denominator term as follows + 1). This expression can be further broken down by the binomial expansion of the term in the brackets into Similarly, the term in the numerator can be binomially expanded to Substituting the expressions for the numerator and the denominator in the expression for H αβ (ḡ, χ) gives One can now notice that (1 − s) −5/2 (1 − t) −5/2 π −3/2 (M α1 M β1 ) −1/2 H αβ (ḡ, χ) would just be a function of (ḡ, χ, s, t) and we can choose to represent it as where A pqrl,1 αβ is a factor that is only purely a function of the masses and temperatures of the species. On comparing the exponents ofḡ 2 and cos χ on both the sides of this expression with these limits, one can also find the summation limits on the the maximum value of r, l in the LHS summation, r ≤ p + q + 1 and l ≤ p + q + 1. Now, on substituting the value of H αβ (ḡ, χ) and H αβ (ḡ, χ = 0) becomes Now transformingḡ into spherical coordinates and integrate over the angles. We also write dΩ = sin χ cos φ and integrate over φ. We also define the Chapman-Cowling cross-section moment integral, in a generalized form, as follows where ζ = k 1/2 The quantity d αβ is an arbitrarily chosen factor that enters the Chapman-Cowling integral, and in general is purely a function of the masses and temperatures of the colliding species, such that d αβ > 0. The coefficients in Zhdanov et al presented in Appendix C are calculated by choosing d αβ = γ α γ β /(γ α + γ β ) (and consequently Now, the coefficient of s p t q in Π αβ is the bracket integral, given by The coefficients A pqrl,1 αβ can be obtained by comparing the binomial decomposition of H αβ (ḡ, χ) with the one solely expressed as a function of (ḡ, χ, s, t), and we ob- where we have usedr andl to indicate that the (r, l) summation on the LHS is different from that in the RHS. Now setting p = r − r 2 − r 4 + l 2 and q = r 2 + r 3 + l − l 2 , and in doing so selecting the s p t q term which fixes l = p + q + r 4 − r 3 − r and l 2 = p + r 2 + r 4 − r, we can write Comparing the limits on both sides, one can show that instead of infinite sum over r, one only needs to sum until r ≤ p + q.

Derivation of
We have Having calculated S p 3/2 (W 2 β )W β , S q 3/2 (W 2 α )W α in the previous section, it becomes much easier to perform the calculation for S l 3/2 (W 2 α )W α , S n 3/2 (W 2 α )W α . Following the previous calculations, we have the bracket integral being the coefficient of s p t q in the generating function Π αα which can then be expressed as Now, we again define a αβ and a βα in the following form such that H αα (ḡ, χ), with the help of the Eqs. (A20) and (A21), becomes which leaves us with needing to derive an expression for W ′ α .W α . In order to do so, we again transform into the variables V α and V ′ α such that where the odd powers ofX will integrate out in H αα .
which on substituting back into W ′ α .W α and defining and on integrating overX, we find So far, this expression agrees with that of Chapman and Cowling. The remaining steps are very similar to the previous derivation, and hence one can find, similar to Eq. (A32) where the N 's and D's are of the form Now, following the binomial expansion for the denominator as in the previous section, one would obtain The binomial expansion for the numerator term remains the same. Combining them, we find the form for H αα (ḡ, χ), Looking at this equation, again it can be expressed as a function in terms ofḡ 2 and cos χ. Therefore, the most general form of Π αα is going to be of the form and on comparing the powers of s and t in both these expressions and selecting the s p t q term, we find that p = r − r 2 − r 4 + l − l 2 and q = r 2 + r 3 + l 2 + l 3 , and in doing so, we inadvertently set l = p + q − r − r 3 + r 4 − l 3 and l 2 = q − r 2 − r 3 − l 3 , because of which the expression for calculating A pqrl,1 where the sum over r need only be computed such that r ≤ p + q, and correspondingly 0 ≤r ≤ p + q + 1 and 0 ≤l ≤ p + q + 1. The general solution is given by This completes the derivation of the second rank-1 bracket integral.

Derivation of
In order to calculate a higher-rank bracket integral, we first notice that fundamentally, only the terms in the curly brackets in Eq. (A37), and in the falling factorial change. Every other term remains the same. Therefore, in the definition of L αβ (ḡ, χ) (defined analogously to H αβ (ḡ, χ), we have Therefore, we need to calculate only the inner product , which can be expressed as Now, W ′ β .W α is given by Eq. (A25), the square of which would be given by And we write down the expression for the product of W ′2 The term in the square brackets is fundamentally what we require; it will replace the one in the curly bracket in the expansion. Additionally, the falling factorial has a 5/2 instead of a 3/2 because we expand a −(r+7/2) αβ from the denominator instead of a −(r+5/2) αβ , with the final expression for the bracket integral given by where the sum over r needs to be only computed till r ≤ p + q, and correspondingly 0 ≤r ≤ p + q + 2 and 0 ≤l ≤ p + q + 2.

Derivation of
Again, as in the previous section, we notice that fundamentally, only the terms in the curly brackets in Eq. (A61) and the index in the falling factorial change. Every other term remains the same. Therefore, in the definition of L αα (ḡ, χ) (defined analogously to L αβ (χ), we have The product (W ′ α .W α ) 2 And we write down the expression for W ′2 Using identity Eq. (A67), we have Therefore, on integrating overX, L αα (ḡ, χ) becomes Following the same procedure as in the previous section, we find Substituting the values of V ′ α .V α and (V ′ α .V α ) 2 − (1/3)V ′2 α V 2 α in the expression for L αα (ḡ, χ), then expanding exp (k αβ − b αβ )ḡ 2 term and shifting the sum for the factors of (k αβ − b αβ )ḡ 2 , we get where the coefficients are given by The term in the square brackets will replace the one in the curly bracket in the expansion for computing A pqrl,2 αα , accompanied by the falling factorial term having a 5/2 instead of a 3/2 because we expand a −(r+7/2) αβ from the denominator instead of a −(r+5/2) αβ . Thus, the general form of the bracket integral and the coefficients A pqrl,2 αα are now given by where the sum over r needs to be computed for r ≤ p+q, and correspondingly 0 ≤r ≤ p + q + 2 and 0 ≤l ≤ p + q + 2.