M Raghunathan

Y Marandet

H Bufferand

G Ciraolo

Ph Ghendrih

P Tamain

E Serre

Multi-temperature Generalized Zhdanov Closure for Scrape-Off Layer/Edge Applications

The derivation of the multi-temperature generalized Zhdanov closure is provided starting from the most general form of the left hand side of the moment averaged kinetic equation with the Sonine-Hermite polynomial ansatz for an arbitrary number of moments. The process of arriving at the reduced higher-order moment equations, with its assumptions and approximations, is explicitly outlined. The generalized multi-species, multi-temperature coefficients from the authors' previous article are used to compute values of higher order moments such as heat flux in terms of the lower order moments. Transport coefficients and the friction and thermal forces for magnetic confinement fusion relevant cases with the generalized coefficients are compared to the scheme with the single-temperature coefficients previously provided by Zhdanov et al. It is found that the 21N -moment multi-temperature coefficients are adequate for most cases relevant to fusion. Furthermore, the 21N -moment scheme is also tested against the trace approximation to determine the range of validity of the trace approximation with respect to fusion relevant plasmas. Possible refinements to the closure scheme are illustrated as well, in order to account for quantities which might be significant in certain schemes such as the drift approximation.

I. INTRODUCTION

Study of plasma transport in presence of impurities in the SOL/edge of current and next-step fusion devices such as WEST, JET and ITER, remains a key topic in magnetic confinement fusion research. Generally, SOL/edge plasmas are treated in the fluid domain by numerical fluid codes usually coupled to kinetic neutrals. The modelling of the friction and thermal forces in such codes is of paramount importance, especially for impurities, as the presence and transport of impurities in the SOL/edge of tokamaks largely depends on the competition between the friction and thermal forces. Generally, for a simple ion-electron plasma of a fixed charge state, there are low-moment closure schemes already available, such as by Braginskii 1 , which help compute these friction and thermal forces using the Landau collision operator. Such a scheme has also been generalized to a higher number of moments in the past [START_REF] Kaneko | Transport Coefficients of Plasmas in a Magnetic Field[END_REF][START_REF] Kaneko | Electrical Conductivity of a High-Temperature Plasma in a Uniform Magnetic Field[END_REF][START_REF] Kaneko | Thermal Conductivities and Coefficients of Thermal Diffusion of a High-Temperature Plasma in a Uniform Magnetic Field[END_REF] for the case of a single-temperature plasma using the linearized Boltzmann operator, and more recently for a multi-temperature plasma using the Landau operator [START_REF] Ji | Closure and transport theory for highcollisionality electron-ion plasmas[END_REF] . These remarkable works have used a large number of moments to achieve convergence of transport coefficients to an increasing degree. However, firstly, there does not seem to be such a work present for the case of multi-temperature plasmas using the Boltzmann operator. Furthermore, the ion-electron case presents a few opportunities to simplify the calculations that cannot be straightforwardly extended to the case of a plasma with multiple species. The aforementioned scheme also solves the reduced set of equations for stationery solutions of plasmadynamical (such as density, momentum density and energy density) and other higher-order thermodynamic quantities (such as heat-flux and the stress tensor). In fluid simulations however, the plasmadynamical quantities are simulated by a time evolution scheme and are not subject to the aforementioned reductions. Therefore, one needs a solution that seeks a relatively station-ary state for the higher-order thermodynamic quantities in terms of the variable plasmadynamical quantities.

In order to remedy these issues, a closure scheme was suggested for the higher-order moments in the linear transport regime [START_REF] Balescu | Transport processes in plasmas[END_REF] , derived from the linearized Boltzmann equation, by Zhdanov et al [START_REF] Zhdanov | Transport processes in multicomponent plasma[END_REF][START_REF] Zhdanov | Diffusion and heat transfer in a multicomponent completely ionized plasma[END_REF][START_REF] Zhdanov | Pfirsch-Schlüter impurity diffusion[END_REF] . The closure scheme, combined with their specific method of obtaining solutions for the higher-order moments such as heatflux and the stress tensor, assuming collisions at common temperature, has come to be known as the Zhdanov closure. The Zhdanov closure has been implemented in a variety of numerical SOL/edge fluid codes, such as Soledge3x-EIRENE 10 , B2-EIRENE [START_REF] Fichtmüller | Multi-Species Developments in the EDGE2D Code[END_REF] , SOLPS [START_REF] Sytova | Derivation of the friction and thermal force for SOLPS-ITER multicomponent plasma modeling[END_REF][START_REF] Makarov | Equations and improved coefficients for parallel transport in multicomponent collisional plasmas: method and application for tokamak modelling[END_REF] , and EDGE2D [START_REF] Bergmann | Implementation into B2 of a 21-moment description for the parallel transport[END_REF] .

Towards this, in the previous article by the authors [START_REF] Raghunathan | Generalized collisional fluid theory for multi-component, multi-temperature plasma using the linearized Boltzmann collision operator for scrape-off layer/edge applications[END_REF] , we derived collision coefficients from the momentaveraged linearized Boltzmann collision operator for a multi-temperature case. The collision coefficients expressed themselves then as linear combinations of product of terms depending on the mass and temperature ratios, and terms depending on the potential of interaction among the colliding particles. The coefficients were then compared to each other, in addition to a single-temperature set provided by Zhdanov et al [START_REF] Zhdanov | Transport processes in multicomponent plasma[END_REF][START_REF] Zhdanov | Diffusion and heat transfer in a multicomponent completely ionized plasma[END_REF][START_REF] Zhdanov | Pfirsch-Schlüter impurity diffusion[END_REF] , and a range of validity was provided for the the singletemperature coefficients. We also illustrated, using the obtained coefficients, a 13N -moment closure scheme in the linear transport regime, showing some of the assumptions behind such a scheme, while establishing the effects of the different sets of coefficients on computed physical quantities such as viscosity and the augmentation in the friction force. We concluded, on the basis of these, that in case of large temperature differences, any set of multitemperature collision coefficients was more reliable than the single-temperature set. Furthermore, it is generally argued [START_REF] Balescu | Transport processes in plasmas[END_REF][START_REF] Zhdanov | Transport processes in multicomponent plasma[END_REF] , that 13N -moments are not adequate for correctly obtaining the values of the transport coefficients. Thus, the question remained open, in the scope of that article, as to whether going beyond the 21N -moment scheme would help. Furthermore, some assumptions behind applying the linear transport regime to higher-order moments were also not explicitly stated. This article seeks to address these topics.

Furthermore, in the past, there have been some attempts at understanding the forces on impurities considering them to be in the trace approximation [START_REF] Braginskii | Transport Processes in a Plasma[END_REF][START_REF] Balescu | Transport processes in plasmas[END_REF][START_REF] Chapman | Thermal Diffusion in Ionized Gases[END_REF][START_REF] Rutherford | Impurity transport in the Pfirsch-Schlüter regime[END_REF] . Some code packages also used to have some trace limit schemes for treating impurities before implementing the 21N -moment single-temperature Zhdanov closure [START_REF] Bufferand | Near wall plasma simulation using penalization technique with the transport code SolEdge2D-Eirene[END_REF][START_REF] Rozhansky | Momentum balance for impurities in SOLPS transport code[END_REF] . Thus, it is also of interest to study the transport coefficients and the thermal and friction forces for the trace approximation against the Zhdanov closure, and thus provide a range of validity for the trace approximation scheme.

In this article, in Sec. II, we introduce the Sonine-Hermite polynomials used, and first re-derive the most general moment-averaged fluid left hand side (LHS) term for an arbitrary number of moments. In Sec. III, we illustrate, step-by-step, the various approximations that reduce the obtained parallel and longitudinal fluid equations to linear forms that can be inverted in order to obtain values for key quantities such as the heat-flux and the stress-tensor, and the friction and thermal forces. We then test the convergence of the scheme with increasing number of moments in IV, comparing the transport coefficients such as the longitudinal viscosity and parallel thermal conductivities, and the convergence of the friction and thermal forces. These tests are performed for a choice of four fusion relevant cases viz. deuteriumtritium as a light multi-species plasma, Carbon, Argon as mid-weight impurities at significant fraction of the main ion species, and Tungsten as heavy impurity in trace quantities. Next, in Sec. V, we consider the special case of the 21N -moment closure and compare it against the trace approximation, establishing a range of validity of the trace approximation for the four chosen cases. In Sec. VI, we seek to illustrate a few extensions to the closure by relaxing some of the assumptions from Sec. III, such as for the linearized Burnett approximation, on the basis of which we build the closure for the drift approximation. The article is finally summarized in Sec. VII. We provide extensive appendices containing the details of the calculations performed, and in particular, provide explicit balance equations in Appendix C, first for the 13N -moments (Appendix C 1), and then balance equations for general vectorial and tensorial moments (Appendix C 2).

II. MOMENT-AVERAGED KINETIC EQUATION

The Boltzmann equation which describes the evolution of the distribution function of a given species α, f α in the frame of the peculiar velocity c α = v α -u, is given by

df α dt + c α .∇f α + 1 m α F * α .∇ cα f α -c αs ∂f α ∂c αr ∂u r ∂x s = β J αβ , (1)
where the common plasma flow velocity u is given by

ρu = α ρ α u α , ρ = α ρ α , (2)
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 force in the moving frame F * α = F α -m α du/dt. The LHS is referred to as the "free-streaming term", and the right hand side (RHS) J αβ is the binary collision term α and every other species of the system. For any quantity ψ α depending purely on species peculiar velocity c α , one can average over Eq. 1 which attains the following form

d dt n α ψ α + n α ψ α ∇.u + ∇.(n α ψ α c α ) -n α dψ α dt + c α .∇ψ α + 1 m α F * α .∇ cα ψ α - c αs ∂ψ α ∂c αr ∂u r ∂x s = R α , (3)
where the ψ α -averaged RHS represents the moment of the collision operator. At this juncture, following previous work, for a multitemperature case, we choose Sonine-Hermite polynomials G mn α of the form,

G mn α (c α , γ α ) = (-1) n n!m α γ -(n+m/2) α × S n m+1/2 γ α 2 c 2 α P (m) (γ 1/2 α c α ), (4)
where γ α = m α /kT α , with the Sonine polynomial S n m+1/2 given by

S n m+1/2 γ α 2 c 2 α = n p=0 - γ α 2 c 2 α p
(m + n + 1/2)! p!(n -p)!(m + p + 1/2)! ,

(5) and the irreducible tensorial monomial p (m) given by

P (m) (γ 1/2 α c α ) = [m/2] l=0 (-2) l γ m/2 α (2m -2l)!m! (2m)!(m -l)! c m-2l α δ l c 2l α , (6)
where the product c m-2l α δ l is symmetrized (see Appendix B for a detailed explanation), such that P (m) is symmetric and traceless for m > 1. The moments b mn α are calculated as

n α b mn α = G mn α f α dc α . (7)
We now substitute ψ α = G mn α in Eq. ( 3), from which we get the the moment-averaged kinetic equation in the form

The collisional RHS R mn αβ is determined from a linearized collision operator using the following ansatz for the distribution function

f α = f (0) α mmax m lm n 2 2n m -2 α γ 2n+m α τ mn b mn α G mn α , (9)
where the distribution function at thermodynamic equilibrium f (0)

α is chosen at the species specific temperature T α , given by

f (0) α = n α γ α 2π 3/2 exp - γ α 2 c 2 α , (10)
and where the constant term τ mn is given by

τ mn = (2m + 1)!(m + n)! n!(m!) 2 (2m + 2n + 1)! .
The number of moments n can be considered to depend on rank-m such that n ≡ n(m). Thus, the number of moments desired can be chosen for each rank-m by fixing the value of l m . E.g. m max = 1, l 0 = 1, l 1 = 0 will lead to the 5N -moment system of equations, m max = 2, l 0 = 1, l 1 = 1, l 2 = 0 will lead to the 13N -moment system, and m max = 2, l 0 = 1, l 1 = 2, l 2 = 1 leads to the 21N -moment system of equations and so on. We shall use the notation n(m) to clarify the rank dependence of the number of moments wherever necessary. In our case following Zhdanov et al, we used the linearized Boltzmann collision operator as follows

R mn αβ ≈ f (0) α f (0) β (G mn α -G mn α )(1 + Φ α + Φ β ) × gσ αβ (g, χ)dΩdc α dc 1β , (11)
where for species α, Φ α is given by

Φ α = (m,n(m)) =(0,0) 2 2n m -2 α γ 2n+m α τ mn b mn α G mn α (12)
This leads to a moment-averaged collision term of the form

R mn αβ = l R mnl αβ , (13)
where are dependent only on the masses, temperatures, number densities, and the potential of interaction between the two species α and β. The expressions for them, alongwith the method of derivation, can be found in the previous article by the authors. In addition, Zhdanov et al provide a set of single-temperature coefficients at the plasma common temperature in Ref. 7-9 for the 21N -moment scheme, which in our previous article, we were able to generalize the calculation for any arbitrary number of moments (we also choose the value of the factor d αβ = γ αβ /2. See Ref. [START_REF] Raghunathan | Generalized collisional fluid theory for multi-component, multi-temperature plasma using the linearized Boltzmann collision operator for scrape-off layer/edge applications[END_REF] and Appendix D for details). For both the multi-temperature and single temperature collision coefficients, we choose the cross section expressions for the shielded Coulomb potential. Note however that the multi-temperature coefficients are generally valid for any temperature differences, but because of the form of the ansatz Eq. ( 9), there in an assumption that the flow velocities of all species are close to plasma common flow velocity, i.e.

R mnl αβ = (1 -δ m0 δ l0 )(A mnl αβ b ml α + B mnl αβ b ml β ) + δ m0 δ l0 C mnl αβ . (14
|u α -u| kT α m α 1/2 . ( 15
)
.

Returning to the LHS, the force is decomposed into general non-velocity dependent body force X, and the electromagnetic force as follows

F * α = X α + Z α e(E + v α × B) -m du dt , (16)
and therefore, the force dependent term can be written as

1 m α F * α .∇ cα G mn α = X αl m α + Z α eE αl m α + Z α e m α {u × B} l - du l dt ∂G mn α ∂c αl , (17)
where the term E + u × B can be thought of as the usual electric field in the moving frame E * . The magnetic field term can be written as

Z α e m α c α × B.∇ cα G mn α = rst ω αt c αs ∂G mn α ∂c αr ( 18
)
where rst is the Levi-Civita tensor and ω α = Z α eB/m α . One can also expand the average of the time derivative term as follows

dG mn α dt = - γ α T α ∂G mn α ∂γ α ∂T α ∂t - γ α T α c α ∂G mn α ∂γ α •∇T α . (19
) It is clear from these expressions that one needs to obtain the various G mn derivatives, contractions, etc beforehand in order to fully expand the LHS. One can use the fact that the irreducible representation of Hermite polynomials is a product of the Sonine polynomial and the irreducible tensorial monomial, and use their properties, such as recurrence, contraction with a vector and rank-2 tensor (1-fold and 2-fold inner products), single and double derivative, to derive the properties for G mn . These properties are reproduced in the appendix A (with some changes for consistency) for the ease of reference.

On using these values and substituting them in Eq. ( 8), the following expression for the balance equation for the general (mn) th moment is found

d dt (n α b mn α )+n α b mn α ∇•u+ ∂ ∂x r n α b m+1,n αr + n γ α n α b m+1,n-1 αr + 2 2m + 1 2m + 2n + 1 2γ α {n α b m-1,n α δ} r + {n α b m-1,n+1 α δ} r - X αr m α + Z α eE r m α + {u × ω α } r - du r dt nn α b m+1,n-1 αr + 2m + 2n + 1 2m + 1 {n α b m-1,n α δ} r + n α ∂u r ∂x s -rst ω αt nb m+2,n-1 αrs + n(n -1) γ α b m+2,n-2 αrs + 2n 2m + 3 δb mn α + n(2m + 2n + 1) (2m + 3)γ α δb m,n-1 α + {δb mn αs } r - 2n 2m -1 δb mn αrs - n(2m + 2n + 1) (2m -1)γ α δb m,n-1 αrs + 2m + 2n + 1 4m 2 -1 2{b m-2,n+1 α δδ} rs + 2m + 2n -1 γ α {b m-2,n α δδ} rs + n α k 2m α n(2m + 2n + 1)b m,n-1 α dT α dt + n α k 2m α n(2m + 2n + 1) b m+1,n-1 αr + n -1 γ α b m+1,n-2 αr + 2 2m + 1 2m + 2n -1 2γ α {b m-1,n-1 α δ} r + {b m-1,n α δ} r ∂T α ∂x r = R mn α = β l [(1 -δ m0 δ l0 )(A mnl αβ b ml α + B mnl αβ b ml β ) + δ m0 δ l0 C mnl αβ ]. ( 20
)
This equation represents the most general fluid moment obtained for a moment average with G mn α . The first line of the equation contains the convective time derivatives d/dt and the space derivatives of the moments (through gradients of scalars and divergences of tensors). The second line contains the velocity-independent force dependent terms. The third and fourth lines contain the viscous-stress and magnetic field dependent terms. The fifth and sixth lines contain the temperature time gradient dT α /dt and space gradient ∇T α dependent terms. The RHS contains the moment averaged collision operator. The general expressions for the RHS R mn αβ for the linearized moment-averaged Boltzmann collision operator for a multi-temperature case, for up to rank-2 moments, can be found in the previous work by the authors [START_REF] Raghunathan | Generalized collisional fluid theory for multi-component, multi-temperature plasma using the linearized Boltzmann collision operator for scrape-off layer/edge applications[END_REF] , which takes the form of a linear combination of moments of a similar rank-m.

In this expression, all terms are symmetrized (See Appendix B), and repeated indices are summed over. It has to be also mentioned that any r th or (rs) th component of a tensor mentioned refers to the additional ranks the tensor has, i.e. the rank-0 quantities do not possess an r th component and rank-1 quantities do not possess an additional (rs) th component. Should such a case arise, the term may be safely set to zero.

The above expression also differs from the general expression Eq. (A1.7) given in Zhdanov et al in Ref. 7, who have a much simpler magnetic field term of the form n α b mn αs rst ω αt δ r . This is just a cosmetic difference, because the Levi-Civita tensor is antisymmetric, and on being contracted with any symmetric term, would vanish. The only term that survives on the expansion is n α rst ω αt {b mn αs δ} r , which can then be expanded out as n α δ r b mn αs rst ω αt . However, we also also have an additional term u × ω term (which could alternately be written as a electric field in the moving frame). Note also the difference in the coefficient of δb mn αrs , and the presence of the Boltzmann constant k multiplying the coefficients of dT α /dt and ∇T α . This expression is similar to the one found in Eq. (3.1.3) of Ref. 20, but defined at the species specific temperature using the Sonine-Hermite polynomials instead of the spherical harmonics.

III. MULTI-TEMPERATURE GENERALIZED ZHDANOV CLOSURE

As one can notice from Eq. ( 20), the moment equation for a general moment b mn of rank-m , contains quantities which are of ranks m ± 1 and m ± 2. This implies that if one truncates the series by choosing a certain number of moments M in the distribution function, one will obtain N corresponding balance equations, but the number of variables in these balance equations will exceed M . This leaves the set of M equations unclosed, requiring elimination of the excess variables through some means. Furthermore, generally most SOL/edge fluid packages solve only for the plasmadynamical quantities, i.e. density, temperature/pressure/energy, and flow momentum density, i.e. (ρ, p/T /E, ρw), which constitute 5N -moments, where N is the number of species. Generally, all additional variables other than these three must be eliminated.

One such direct method is the closure recommended by Grad himself, which involves calculating the higher order moments in terms of lower ones by using the ansatz for the distribution function Eq. ( 9) in the expression for the moment Eq. (7). For symmetric and irreducible moments, assuming the ansatz of the distribution function does not contain the moment, such a process is equivalent to setting to zero that moment. It has come to be known as Grad's closure in literature, and is the finishing touch in the description of Grad's method. This generally works well for fluids such as monoatomic gases, where the RHS vanishes in the hydrodynamic equations and only the LHS needs to be treated. In particular, in a 5N -moment approximation, the heat-flux moment b 10 α , i.e. h α = 0 (thus the conventional heat flux q α = (5/2)p α w α ) and the stress tensor π α = 0 from Grad's closure. One can notice that it does not recover the temperature-gradient force term in RHS of the momentum balance equation, because heat-flux b 11 α calculated in this manner does not depend on the temperature gradients. It also does not recover the usual visco-elastic form of the stress tensor, where the stress tensor is proportional to the rate-of-strain tensor. Thus, this leaves us with an oversimplified description of the heat-flux and the stress tensor if the 13N -moment Grad's closure is used. Also notice that this closure typically leads to the higher-moment of any species depending only on the lower moments of that species only.

One other method was proposed by a series of authors [START_REF] Braginskii | Transport Processes in a Plasma[END_REF][START_REF] Kaneko | Transport Coefficients of Plasmas in a Magnetic Field[END_REF][START_REF] Kaneko | Electrical Conductivity of a High-Temperature Plasma in a Uniform Magnetic Field[END_REF][START_REF] Kaneko | Thermal Conductivities and Coefficients of Thermal Diffusion of a High-Temperature Plasma in a Uniform Magnetic Field[END_REF][START_REF] Balescu | Transport processes in plasmas[END_REF][START_REF] Landshoff | Convergence of the Chapman-Enskog Method for a Completely Ionized Gas[END_REF][START_REF] Spitzer | Transport Phenomena in a Completely Ionized Gas[END_REF] , under what is known as the linear transport assumptions, to find approximate equations for the higher-order moments, usually using the first Chapman-Enskog approximation, and ignoring some additional terms on the basis of low electron mass, for a single-temperature plasma, to find transport coefficients for an ion-electron plasma. These usually involve a simultaneous solution of the approximate linearized equations for ions and electrons, which after the approximations may or may not be coupled, to find the transport coefficients. Some of these works use the linearized Landau collision operator, and the others the linearized Boltzmann operator. These works furthermore find the usual form of the friction and thermal forces, in terms of the flows and the temperature gradients respectively. However, the transport coefficients and the forces calculated from such a scheme, which is only strictly applicable to an ion-electron single-temperature plasma, are not usually applicable to a plasma with impurities at different charge states in significant amounts. Thus, to address this problem, Zhdanov et al proposed a new solution scheme for the case of impurities in a plasma.

A. General assumptions

Zhdanov et al proposed a new scheme [START_REF] Zhdanov | Diffusion and heat transfer in a multicomponent completely ionized plasma[END_REF][START_REF] Zhdanov | Pfirsch-Schlüter impurity diffusion[END_REF] that involves linearizing the balance equation for the higher-order moments under certain assumptions as follows 1. Firstly, a desired set of moments is chosen, comprising of plasma dynamical moments (ρ, ρw, nkT ), thermodynamically privileged higher-order moments (h, π), and thermodynamically non-privileged moments of even higher-order. The balance equations for this set of moments is calculated with the collision terms being calculated by the ansatz of the distribution function containing these moments Eq. (9).

2. For moments out of this desired set, Grad's closure is used on them, i.e. for our case of symmetric irreducible moments, these higher-order moments are set to zero.

3. The plasma dynamical moment (ρ, ρw, nkT ), are considered zeroth order in K n . Their space gradients are considered again to go as first order in K n , where

K n ∼ λ L sc , τ τ sc , (21)
where λ represents mean free path between collisions, and τ represents mean time between collisions, and L sc and τ sc represents the scale lengths and the scale times of the system in question. However (w, T ) are higher than first order, the products of (w, T ) with themselves or other moments may be considered to be of order higher than one in K n .

4. The higher-order moments, both privileged and non-privileged, are considered to be of the order of one or higher in Knudsen number K n .

5. The time derivatives of higher-order moments, both privileged and non-privileged, are neglected, meaning that quantity that the higher-order moment represents changes slowly over the characteristic timescale of the system τ sc . This essentially means that the time derivatives of higherorder moments are considered larger than first order in Knudsen number K n . Physically, it means that the moment evolves slower than any changes in the moment caused by collisions.

6. The space gradients of non-privileged higher-order moments are also neglected, which means that the macroscopic quantities represented by the higherorder moments change gently over the scale length of the system L sc , meaning the gradients of the higher-order quantities are of an order larger than one in K n .

7. Products of a higher-order moment with other moments, higher or lower order, are also neglected, since the product of moments which are the first order in Knudsen number will lead to quantities which are second order or higher in K n .

8. The common flow u is considered zeroth order in K n , but its space and time gradients are considered between order zero and one in K n . 9. Consequently, γ can be considered to go as order one, and factors of γ -1 may be considered to be between (-1) st in K n . Factor of k/m may be considered even lower than the (-1) st order. (Which would make some gradients like γ -1 ∇w and (k/m)∇T equal to or lower than order one.)

These approximations essentially transform the balance equations of these higher-order moments into linear, nondifferential equations, which can then be solved to obtain approximate values of these higher-order moments. It is similar to the first-order approximation of the Chapman-Enskog scheme, because of the explicit ordering in terms of the Knudsen number. (Note however that it is not exactly the same as the Burnett approximation. We shall touch lightly on this point later). A similar scheme has also been proposed by Balescu 6 , which retains the time derivatives of the higher-order moments however, in addition to the explicit ordering in terms of a provisionally defined hydrodynamic timescale τ H rather than the Knudsen number.

B. Specific assumptions in the reduction of the balance equations

On applying the approximations as outlined, we first find that the balance equations for the plasmadynamical quantities (ρ, ρw, nkT ) survive as it is, as can be seen from Appendix C 1. However, in the LHS of these balance equations, we find that there are open variables such as the heat-flux h α and the stress tensor π α . Again, on applying the approximations from the previous subsection, we find the balance equations for the heat-flux h α and the stress tensor π α reduce to the following

-π αl X αl m α + Z α eE αl m α + {u × ω α } l - du l dt - 5 2
kT α m α ∇.π α + 5 2 k m α n α kT α ∂T α ∂x r -h α × ω α = β R 11 αβ , (22)
and

4 5 5 2 n α kT α ∂w αr ∂x s + ∂h αr ∂x s -ρ α w αr X αs m α + Z α eE αs m α + {u × ω α } s - du s dt + 2n α kT α ∂u r ∂x s -2{π αsl rst ω αt } = β R 20 αβ (23)
respectively. Furthermore, from the RHS of the balance equations, we can observe that there are higher-order moments which remain unclosed. Therefore, the balance equations for the general higher-order vectorial and tensorial non-privileged moments b can be similarly reduced to

-n(1)n α b 2,n(1)-1 αl X αl m α + Z α eE αl m α + {u × ω α } l - du l dt - 2n(1) + 3 3 n α b 0,n(1) α X α m α + Z α eE m α + u × ω α - du dt -b 1n(1) α × ω α = β R 1n(1) αβ , (24)
and

- 2n(2) + 5 5 n α b 1,n(2) α X α m α + Z α eE m α + u × ω α - du dt -2{b 2n(2) αsl rst ω αt } = β R 2n(2) αβ (25)
In the absence of body forces X α and parallel electric fields E , and a static common flow, these equations for the heat flux h α and stress tensor π α can be further resolved along the direction parallel to the magnetic field, such that we obtain

- 5 2
kT α m α ∇ .π α + 5 2 k m α n α kT α ∇ T α = β R 11 αβ , (26)
respectively, where for a vector a, a = bb.a and for a traceless symmetric tensor A, A = (bb -δ/3)(bbδ/3) : A, where b is a unit vector along the direction of the magnetic field. Now, one can further assume that the parallel gradients of the higher-order moments are weak in nature, and reduce these equations further down to

5 2 k m α n α kT α ∇ T α = β R 11 αβ , (28)
and

2n α kT α ∂u r ∂x s = β R 20 αβ . (29)
(Such an assumption may not be valid under certain conditions, to which we will return to in Sec. VI). Under the same assumptions of no body forces and no parallel electric fields, the balance equations for the non-privileged higher-order moments in the parallel direction can be reduced as follows b

1n

(see Appendix C 2) can be similarly reduced to

0 = β R 1n(1) αβ , (30)
and

0 = β R 2n (2)
αβ .

The expressions Eqs. ( 28)-( 31) represent then a closed set of equations which form the basis of the generalized Zhdanov closure.

C. Solution of the linear reduced system of balance equations

We now proceed to explain how to solve the linear system of equations given by Eqs. ( 28)- (31). First, we define M mnl αγ such that

M mnl αγ = A mnl αα + B mnl αα + β =α A mnl αβ , α = γ B mnl αγ , α = γ ( 32
)
where the index γ runs over all N species. It is also similar to the notations q nl , p nl used in Refs. 23-27. We then define an

N × N matrix M mnl M mnl =       M mnl αγ M mnl αδ . . . M mnl αω M mnl δα . . . . . . . . . . . . . . . M mnl ωα . . . . . . M mnl ωω       , (33)
such that

R mn α = l M mnl α B ml α , (34)
where M mnl α is the row-α, and where B mnl is a column vector of length N given by

B ml =      b ml α b ml γ . . . b ml ω      . ( 35
)
Note that, for the special case of m = 0, l = 1,

R 0n α = β C 0n0 αβ + l=1 M 0nl α B 0l α . (36)
For the rank-1 quantities, we define an (l 1 -1)×(l 1 -1) block matrix such that,

M 1 l1 =       M 111 M 112 . . . M 11l1 M 121 . . . . . . . . . . . . . . . M 1l11 . . . . . . M 1l1l1       , (37)
and a block column vector given by

M 1 0l1 =      M 110 M 120 . . . M 1l10      , (38)
such that the linear transport relations can be represented as

Λ 1 l1 T = M 1 0l1 W + M 1 l1 H 1 l1 , (39)
where H 1 l1 and Λ 1 l1 are block column vectors of length (l 1 -1) given by

H 1 l1 =      H B 12 . . . B 1l1      , Λ 1 l1 =      Λ 0 . . . 0      , ( 40
)
where

T , W , H , B 1l , are column vectors of length N containing values of ∇ T γ , ρ γ w γ , h γ , b 1l γ respectively, and Λ is a N × N diagonal matrix containing 5 2 k mγ n γ kT γ on its diagonal.
Similarly, for rank-2 quantities, a l 2 × l 2 block matrix M 2 l2 can be defined as follows

M 2 l2 =       M 200 M 201 . . . M 20l2 M 210 . . . . . . . . . . . . . . . M 2l20 . . . . . . M 2l2l2       , (41)
from which the set of linear transport relations for rank-2 quantities can now be represented as

-2P 2 l2 = M 2 l2 Π 2 l2 ( 42
)
where Π 2 l2 and P 2 l2 are block column vectors of length

l 2 Π 2 l2 =      Π B 2l
. . .

B 2l2      , P 2 l2 =      P 0 . . . 0      , (43)
where Π , P are column vectors of length N containing values of π γ and n γ kT γ respectively. The symbol is shorthand for {∇u}.

The solutions to Eqs ( 39) and ( 42) are given by inverting the equations as follows

H 1 l1 = (M 1 l1 ) -1 Λ 1 l1 T -(M 1 l1 ) -1 M 1 0l1 W (44)
Π 2 l2 = -2(M 2 l2 ) -1 P 2 l2 . ( 45
)
Thus one can obtain expressions for rank-1 general moments in terms of the flow velocities w and the temperature gradients ∇ T , and expressions for rank-2 such that it is expressed in terms of the longitudinal rate-of-strain tensor 2 . One can confirm that the higher order moments calculated in this manner respect the orderings as mentioned in the previous subsection, for example, from Eq. (8.4.6) of Ref. 7, where the heat flux term h α is a factor of γ -1 α τ α smaller than the flows and the temperature gradients, i.e. one order higher in K n . The gradients of such terms will be one order higher than K n because of the assumptions on the gradients of the plasmadynamical quantities. Thus, the ordering chosen is generally respected by the results obtained from the closure.

Furthermore, the partial parallel thermal conductivities λ δγ can be found in the (1, 1) element of the block matrix (M 1 l1 ) -1 Λ 1 l1 and similarly, the partial longitudinal viscosities η γ can be found in the first element of the block column vector (M 2 l2 ) -1 P 2 l2 . The full parallel thermal conductivities λ δ and full longitudinal viscosity η can be obtained by summing over the columns of the respective matrix elements, i.e. λ δ = γ λ δγ , and η = γ η γ respectively.

A version of this closure for 13N -moment case was illustrated in the previous article by the authors, however, neglecting the self-collisions, for a purely non-magnetic case [START_REF] Raghunathan | Generalized collisional fluid theory for multi-component, multi-temperature plasma using the linearized Boltzmann collision operator for scrape-off layer/edge applications[END_REF] . The non-magnetic case is superficially equivalent to the case here where we only consider parallel and longitudinal transport coefficients.

D. Friction and thermal forces

Once the values of the higher-order rank-1 moments have been calculated, they can be substituted back into the RHS block of the momentum equation, i.e. the balance equation for ρ α w α , for all species as follows

R 00 = M 100 W + M 1 l10 H 1 l1 ( 46
)
where H 1 l1 is given by Eq. ( 44), and a block row vector M 1 l10 is given by

M 1 l10 = M 101 M 102 . . . M 10l1 . (47)
On substituting H 1 l1 , we get

R 00 = [M 100 -M 1 l10 (M 1 l1 ) -1 M 1 0l1 ]W +M 1 l10 (M 1 l1 ) -1 Λ 1 l1 T . ( 48
)
The part of the collisional force dependent on the flows W is termed the "friction force" and the part dependent on the temperature gradients T is termed the "temperature-gradient force" (or, "thermal force" at times) [START_REF] Braginskii | Transport Processes in a Plasma[END_REF][START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF] . As one can observe, the addition of terms from H 1 l1 "augments" the friction force term depending on W , which has more contribution the higher the number of moments chosen.

IV. CONVERGENCE OF THE GENERALIZED ZHDANOV CLOSURE

In order to test the closure scheme, following the previous article, we investigate four cases of three-component plasmas, i.e. electrons and two other species, the most relevant to SOL/edge physics, as follows

• The fusion fuel, containing deuterium and tritium (D-T), at comparable densities, as D-T fusion is planned to be used in current and future burning plasma campaigns,

• Light impurities at significant fraction (10%) of the main fuel species, i.e. hydrogen and carbon (C-H), with the carbon in the plasma originating from facing plasma components made of graphite,

• Injected mid-weight impurities with densities at a small fraction (1%) of the fuel species density, e.g. hydrogen and argon (Ar-H), often used for controlled experimentation with impurities, for measurement of background plasma flows, etc [START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF] , or for other purposes

• Heavy impurity at trace levels (0.001%), i.e. hydrogen and tungsten (W-H), where the tungsten usually originates from the walls and divertors made of tungsten.

The parameters chosen for these cases can be found summarized in Table I.

A. Convergence of transport coefficients

For convergence calculations, we choose the maximum number of rank-1 and rank-2 moments to be l 1 = 8 and l 2 = 8 respectively, i.e. five and six more moments each in addition to r α and σ i respectively. The choice for this is mainly motivated by the compute time in Mathematica [START_REF] Inc | [END_REF] .

First, we graph the total longitudinal viscosity and total parallel heat conductivities for the D-T case against the temperature ratio for the multi-temperature coefficients with d αβ = γ αβ /2, and they can be found in Fig.

( 1) and ( 2) respectively. One can notice that the total longitudinal viscosity η converges remarkably quickly. The addition of a second tensorial moment, i.e. M = 2, results in a significant jump, but the addition of more moments simply seems to make the curves of the total longitudinal viscosity overlap. Generally, we observe the same for the total parallel thermal conductivities, however, addition of a third moment in addition to r α seems to weakly contribute to the calculation of λ T . This is in agreement with the results found by Balescu 6 (Fig. 4.1, pp. 238), that addition of more moments beyond those in the 21N -moment range seems to make no significant changes in these transport coefficients. Furthermore, the values of the longitudinal viscosity and parallel thermal conductivity values are well within the order of magnitude values as of Braginskii At this point, a question may be posed as to why we are plotting quantities with respect to the temperature ratio and not β α = ω α τ αα as has been done traditionally in Refs. 2, 5, and 6. In these references, a form of the collision operator is used which approximates the cross section such that the collision coefficients are essentially some rational multiple of τ αα , which is proportional to the Coulomb logarithm. Thus, on forming the matrix of collision coefficients and inverting the matrix, one essentially finds transport coefficients to be a polynomial of β α . However, as one can see from our collision cross section, there are other terms in addition to the logarithmic term, some of which may be significant compared to the Coulomb logarithm. This fact does not allow us to claim that our transport coefficients are merely polynomials of β α , making it an inconvenient plotting parameter. Furthermore, since we are only looking at the parallel/longitudinal component, the strength of the magnetic field is inconsequential to these parameters. We, therefore, choose the temperature ratio along with the distinction of the four physical cases, hoping to shed light on them individually. However, generalization of the multicomponent parallel/longitudinal closure to magnetized closures will comprise a part of our future work, and we may address this problem again later in that context.

In order to study these convergence effects more quantitatively, we calculate the differences from references values. The first is the maximum difference of the transport coefficients obtained with the multi-temperature coefficients with respect to the curve obtained for the same multi-temperature moments with maximum number of moments. We call this difference I, which allows us to see how quickly the transport coefficients from the multi-temperature scheme converge. The second difference II, is the maximum difference of the transport coefficients obtained with single-temperature coefficients with respect to the values obtained with the maximum number of moments using the single-temperature coefficients, i.e. M = 7 for rank-1 coefficients and M = 8 for rank-2 coefficients. This allows us to see how the transport coefficients in the single-temperature scheme converge. The third difference III we define is the difference of the single-temperature transport coefficients with respect to the values for the transport coefficients obtained from the multi-temperature scheme for the same number of moments, which allows us a global perspective of differences between the single and multi-temperature schemes. Now we calculate these difference values for the different transport coefficients using the multi-temperature and single-temperature coefficients for the four physical cases chosen, which can be found in Table II. Again, to avoid difference values being exaggerated by transport coefficients approaching zero, we restrict the temperature ratio range to 0.5-2.

From the table, a few observations can be made 1. The most significant difference seems to indeed remain between the 13N -moment scheme and the others. An improvement in the range of 10 -150% can be seen, with the lighter impurities representing the lower end of the range, and the mid-weight impurities representing the higher end.

2. The total longitudinal viscosities seem to converge much faster than the total thermal conductivities in general. All viscosities seem to converge to two decimal places in just the second moment (i.e. the 21N -moment scheme).

3. In general, the differences for any transport coefficient seem to increase as the weight of the impurity increases and attain a maximum between the Carbon and Argon cases.

4. The thermal conductivity λ β (β = T, C, Ar, W ), seems to not even converge to one decimal place for Carbon and Argon within the 21N -moment scheme.

The thermal conductivity λ β (β = T, C, Ar, W ), seems to converge better than λ α , but still exhibits significant differences for Carbon and Argon. Additionally, one can roughly say that the thermal conductivities converge to one decimal place in the 21N -moment scheme are not converged to one decimal place for mid-weight impurities.

5. The thermal conductivity λ e seems to converge to two decimal places in the 21N -moment scheme even for mid-weight impurities, in contrast with the other two thermal conductivities λ α and λ β . Furthermore, it seems to converge to two decimal places in the scope of the 21N -moment scheme, slightly better than the latter conductivities. This seems to corroborate why extension of the 13Nmoment scheme to the 21N -moment scheme was considered sufficient for the λ e in previous works.

6. Differences of types I and II seem to more or less be of the same order for any considered moment.

7. From the type III columns, we can observe that the differences between the single-temperature scheme and the multi-temperature scheme remain quite significant at any number of moments, and remain virtually the same beyond the 21N -moment scheme. The average difference between these two schemes seems to be in the range of 10 -80%, with the mid-weight impurities representing the higher side of this range. Furthermore, even though it might seem as if they type II errors are converging marginally faster for the mid-weight impurities, significant differences still remain compared to the multi-temperature scheme (type III), even for heavy impurities. Based on this, at least for midweight impurities and heavy impurities, we recommend using multi-temperature coefficients.

B. Convergence of the friction and thermal forces

A corollary, but important aspect to the effect of addition of more moments on the transport coefficients, is its effect on the finally calculated values of the friction and temperature-gradient dependent forces in the RHS of the momentum equations. This is especially so, because the fluid code packages mentioned in the previous section explicitly use the Zhdanov closure scheme to calculate their friction and thermal forces, making it necessary to study the implementation of the scheme. The augmentation of the friction force can be studied, for any moment with respect to the 5N -moment scheme, by the following expression

% f ric,aug (l 1 ) = M 1 l10 (M 1 l1 ) -1 M 1 0l1 M 100 × 100%, (49)
where, since the M -matrices are N ×N in dimension, the division is performed element wise. Here we use the term "augmentation" in a general sense -additional moments may also contribute to reducing the values of the forces. We studied the augmentation of the friction force for the case of l 1 = 1 in the previous article [START_REF] Raghunathan | Generalized collisional fluid theory for multi-component, multi-temperature plasma using the linearized Boltzmann collision operator for scrape-off layer/edge applications[END_REF] . The addition of the thermal force can be defined with respect to the first contribution from the heat-flux terms., i.e. w.r.t. the 13N -moment scheme (since there is no temperature gradient force in the 5N -moment scheme), as follows

% therm,aug (l 1 ) = M 1 l10 (M 1 l1 ) -1 Λ 1 l1 -M 1 10 (M 1 1 ) -1 Λ 1 1 M 1 10 (M 1 1 ) -1 Λ 1 1
×100%.

(50) One element in these matrices of special interest, is the effect of the background flow on the impurity, which in this case is the effect of deuterium flow on the tritium friction force term, since this term is what affects the impurity dynamics.

We first proceed to plot, as in the previous sections, the augmentations of the friction and thermal forces for the deuterium-tritium case for this element, in Figs. 3 and4 respectively. One can see that the forces are converged relatively quickly on the addition of just the second moment, i.e. under the 21N -moment scheme. However, a key difference between the singletemperature and multi-temperature scheme is that the single-temperature scheme suggests very little variation in friction forces over the temperature range, whereas the multi-temperature scheme suggests a much larger range of magnitude of friction forces over the chosen temperature range. A similar observation can be made for the thermal forces as well. Generally, the friction and thermal forces for the deuterium-tritium case seem converged adequately near the unity temperature ratio for the 21Nmoment case, but diverge as the temperature ratio goes further away from one.

Next, as in the previous subsection, we form tables of augmentation of the friction and thermal forces as defined in Eqs. ( 49) and ( 50) for the four different physical cases chosen. As before, we define three types of augmentations. Type I refers to the augmentations for the multi-temperature case, Type II represents augmentations for the single-temperature case, and Type III denotes the percentage difference of the single-temperature force w.r.t. the multi-temperature force for any given moment. As mentioned earlier, we consider the force coefficients for the forces on the impurity due to the main-ion flow velocity. The tables for the friction force augmentation and the the thermal force augmentation can be found in Tables III and IV. From the table of the friction forces, we can generally observe that irrespective of the single-temperature or multi-temperature coefficients, the friction force seems to converge to two decimal places within the 21N -moment approximation compared to the friction force of the 5Nscheme. In fact, the 13N -moment scheme seems to increase the friction force most significantly, with the 21Nmoment scheme delivering further precision. We can also observe that there still remain significant differences between the single-temperature and multi-temperature friction forces at any given number of moments, which does not decrease significantly on addition of more moments. They also seem to remain more-or-less of the same order for increasing impurity weight. However, a part of this effect is counteracted by the fact that the augmentations themselves seem to be very small for heavy impurity case of tungsten, about 1% or less. However, for mid-weight impurities such as carbon and argon, using multi-temperature coefficients may be necessary for more precise representation of the friction force, as the augmentations to the friction force remain significant.

In mild contrast, the thermal forces generally seem to converge to one decimal place or less within the 21Nmoment scheme, however, the addition by the 21Nmoment scheme is significant for low-weight and midweight impurities. For the deuterium-tritium and carbon case, the 21N -moment thermal force augmentation seems fairly close to the values on addition of further moments. However, for argon and tungsten, the convergence of the augmentations of friction and thermal forces is much slower. For tungsten however, the friction and thermal forces themselves remain much smaller compared to mid-weight impurities in higher concentrations, as can be noticed from Fig. 5, therefore, a slow convergence does not affect the transport coefficients related to it significantly. Furthermore, similar differences between the multi-temperature and single temperature coefficients, as seen from the Type III difference, remain significant. (The ones for the D-T case are slightly exaggerated by thermal force being close to zero for a certain temperature ratio. But from Fig. (3), one can see that the differences are still significant).

On the basis of the studies on the friction force and the thermal force, we find generally that the 21N -moment scheme brings added precision to the force values over the 13N -moment scheme. Again however, we advise caution while using the single-temperature coefficients for midweight impurities, since the convergence of the forces is much lesser than that of low-weight impurities, and since the force augmentation is much more significant than heavy impurities.

C. Section summary

Based on the observations in the two previous subsections, we can generally say that the 21N -moment system represents a significant gain in precision of transport coefficients over the 13N -system. More precisely, the 21Nmoment system can represent the total viscosity very well (to two decimal places), but the total thermal conductivities are represented only moderately well (to maximum one decimal place). Furthermore, there still persist significant differences between the single-temperature coefficients and multi-temperature coefficients, particularly for mid-weight impurities and heavy impurities. Consequently, even the converged values of transport coefficients for multi-temperature and single-temperature coefficients differ significantly. Similar observations are made for the convergence of the friction forces and the thermal forces are made as well, with the friction force being precise to two decimal places in the 21N -moment scheme, and the thermal forces to one decimal place, with similar reservations about the convergence for midweight impurities. Thus, caution is recommended when using the 21N -moment single-temperature scheme for mid-weight impurities such as Carbon and Argon, and heavy impurities such as tungsten. For them, the multitemperature scheme may be better suited. We have thus, with the help of the convergence study, managed to establish the significant improvement brought over by the multi-temperature 21N -moment scheme, and have also managed to establish a range of validity for the existing 21N -moment single-temperature which comprise the Zhdanov closure.

V. THE 21N -MOMENT ZHDANOV CLOSURE SCHEME VS TRACE APPROXIMATION

At this juncture, it would be interesting to study the performance of trace approximation schemes with respect to 21N -moment Zhdanov closure, because the earlier versions of some of the fluid code packages often depended on some sort of trace approximation on the impurities for simulations. Having demonstrated in the previous section that the 21N -moment schemes bring the values of transport coefficients, and, friction and thermal forces to two decimal places of convergence, we can use it as a standard to study the accuracy of any trace approximation scheme. The advantage of using trace approximation in the Zhdanov closure scheme is that they are computationally lighter as compared to the inclusion of all collisional coefficients. This is an important consideration, as generally in the numerical implementation of the Zhdanov closure, the calculation of the collision coefficients is what takes the most amount of time, the inversion of the matrices, e.g. in Eqs. ( 44) and ( 45), being relatively quick [START_REF] Bufferand | Three-dimensional modelling of edge multicomponent plasma taking into account realistic wall geometry[END_REF] . Thus, if the domain of validity of the trace approximation is established, one could identify certain scenarios in which the fluid code could switch to the trace approximation to save compute time.

In order to apply the trace approximation in the Zhdanov closure, the impurities are assumed to be at such low concentrations that they have no effect on the background plasma at work. Thus, one way to define the trace approximation as follows 1. The self-collisions for impurities are excluded by the means of excluding the (A mnl αα + B mnl αα ) term in Eq. ( 32) in M mnl αβ .

2. Terms A mnl αβ , B mnl αβ are neglected for electronimpurity collisions, i.e. we neglect any effects collisions with impurities have on electrons.

3. Furthermore, we neglect A mnl αβ , B mnl αβ for the ionimpurity collisions, i.e. we assume the impurities are not in a sufficient quantity to affect the main ion species.

Then this trace approximation assumptions are used for the collision coefficients which enter the 21N -moment Zhdanov closure, and the transport coefficients and the friction/thermal forces are calculated in that manner. In principle, about 33% of compute time can be saved by ignoring these three coefficients.

On comparing the transport coefficients and the friction and thermal forces for the multi-temperature trace tungsten case, we find that the deviations (excepting for λ W , which generally has a negligible contribution to the heat-flux for trace density values) are less than 1% and 0.1% respectively. This is expected, as such trace approximations are usually only applicable to impurities at very low concentrations, such as tungsten being present at 10 -5 the density of the main ion species. They are generally valid in either fluid codes where there is some uncertainty of values of transport coefficients between impurities and the main species, when on doing the trace approximation at least would ensure no spurious contributions from the impurity to the main plasma. They are also used in some orbit-following codes where it is hard to quantify the back-reaction of the impurities on the plasma, because the main plasma is taken to be effectively static. Generally, such codes are unable to simulate impurities at significant concentrations.

Thus, it is worth studying the transport coefficients for the three other physical cases mentioned relevant to SOL/edge plasmas as in the previous section. As in the previous section, we calculate and tabulate the differences from a reference curve. We can safely take the the 21N -moment, multi-temperature values with full collisions included, for each physical case, as the reference, and define differences with respect to this curve, defined as follows

% dif f, case = t case -t 21N,multi-temp,f ull t 21N,multi-temp,f ull × 100%.
(51) The of the maximum percentage errors for different cases can be found in Table V.

From this table, we can find a few general patterns. Firstly, the differences in transport coefficients seem to increase initially with impurity weight/charge and then decrease. This is consistent with our observations in our previous article [START_REF] Raghunathan | Generalized collisional fluid theory for multi-component, multi-temperature plasma using the linearized Boltzmann collision operator for scrape-off layer/edge applications[END_REF] , that the differences attain a maxima in between the charges/densities of Carbon and Argon. This is partly consistent with increasing differences observed by Balescu et al 6 with increasing charge state Z α , who did not observe a decrease because they did not consider the densities to decrease as the charge state increased. Thus, the decrease observed in our case arises mainly from the decreasing densities across our four physical cases. And second, the transport coefficients calculated from single-temperature coefficients exhibit larger differences from those calculated from their multi-temperature counterparts.

Concerning the friction and thermal forces, as in the previous section, we compare the coefficient of the force on the impurity due to the background ion flow or temperature gradient. However, for low and mid-weight impurities, we find that though friction forces for all cases lie within 30% of each other, the thermal forces are overestimated by the trace approximation by a factor of two. The maximum difference is observed for the carbon case, similar to what was observed earlier for the case of transport coefficients, as can be seen in Tables VI andVII. We can also notice that the single-temperature coefficients present significant differences compared to the multitemperature ones, in agreement with all our observations so far.

In order to illustrate better the differences of using the trace approximation in a case which does not allow for it, we proceed to plot the values of the transport coefficients for multi-temperature and single-temperature cases, for the full collisions and trace approximation cases. Fig. 6 represents the plot of the total longitudinal viscosity and total parallel thermal conductivities, for the trace approximation and full closure, for the D-T case. From the plot for viscosity and thermal conductivities, one can indeed notice that the trace approximation overestimates the transport coefficients significantly. Furthermore, the curves with the single-temperature and multitemperature schemes still follow each other quite closely, only separating out significantly for trace approximation.

The thermal conductivities seem to be closer to each other, with a difference of about 30 -50% between the trace approximation and full collisions, but the viscosity is overestimated by a factor of two nearly. One can also notice from Fig. 7, that indeed though the friction forces follow each other to some 30% difference, the thermal force is overestimated by a factor of two nearly.

It is worth mentioning that there is also a heuristic manner in which the trace approximation can be formulated, which lends itself to a much faster computational scheme than the abovementioned full and trace Zhdanov closure. This involves the same assumptions on the collision coefficients, but then proceeds to formulate the friction and thermal forces on the particles by the means of solving the thermal diffusion of impurities from the momentum equation at slow time evolution and small spatial gradients approximation [START_REF] Chapman | Thermal Diffusion in Ionized Gases[END_REF][START_REF] Rutherford | Impurity transport in the Pfirsch-Schlüter regime[END_REF][START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF] , with the collision coefficients at common plasma temperature. Generally, this procedure gives a heuristic estimate on the coefficients of the temperature gradients in the thermal force, where both coefficients are generally proportional to the square of the impurity charge state Z 2 α , and for heavy impurities, the ion thermal gradient coefficient is of the same order to the electron one, such that their ratio is of order O(1). If the thermal force is such that

F therm ∼ α e ∇ T e + β i ∇ T i , (52)
then the values of α i and β e are given by

α e = α e0 Z 2 (53)
β i = 3 µ + 5 √ 2Z 2 (1.1µ 5/2 -0.35µ 3/2 ) -1 2.6 -2µ + 5.4µ 2 (54)
where µ = m Z /(m i + m Z ) the relative mass ratio of the impurity to the main ion, and Z is the charge state of the impurity. This heuristic scheme can also be considered as a formulation of the thermal force in the static 5Nmoment single-temperature approximation, which is also the fastest in terms of performance. According to Ref. 17 and 28, α e0 = 0.71, and according to Ref. 16, α e = 1.0. Thus, a range of heuristic thermal force scaling factors β i /α e can be defined of the form

(3, 4.23) × µ + 5 √ 2Z 2 (1.1µ 5/2 -0.35µ 3/2 ) -1 (2.6 -2µ + 5.4µ 2 )Z 2 , ( 55
)
where this scaling factor measures the relative magnitude of the forces from the ion and electron temperature gradients. We can state that for the tungsten case at unity temperature ratio, the scaling is followed very well, as is expected at the trace limit. But there exist significant deviations outside the vicinity of equal temperatures. Furthermore, there are deviations for mid-weight impurities, where this heuristic approximation overestimates the scaling compared to Zhdanov closure in the trace limit. Therefore, such a heuristic thermal force scheme is not recommended for mid-weight impurities at significant densities. This scheme is only strictly valid for the trace impurity case, and where the temperatures of all species are same, which it were to be the case would be the least computationally expensive.

Thus, from the observations in this section, we can state confidently that the trace approximation by neglecting certain collisional coefficients is indeed valid for the trace heavy impurity case, while exhibiting significant differences in the transport coefficients and the thermal force for the mid-weight impurities at significant density fraction. Thus, on this basis, we can recommend safely replacing the 21N -moment multi-temperature closure by its trace counterpart for heavy impurities at the trace limit, like for the tungsten case. Furthermore, one can also use the heuristic scheme for the thermal forces if the temperature ratio is close to unity. At this point, we mention a caveat that all treatment so far is only valid when radiative energy loss can be neglected, which may not be the case for heavy impurities even at low concentrations. Thus care must be taken to not use the any of the trace or full schemes in cases where radiation losses are significant. Now that we have tested the 21N -moment multitemperature Zhdanov closure for its convergence and performance against the trace approximation, we now proceed to relax some of the assumptions in order to illustrate a few specific cases of application.

VI. RELAXING SOME ASSUMPTIONS ON THE 21N -MOMENT MULTI-TEMPERATURE CLOSURE

The assumptions mentioned in Sec. III generally are valid for the linear transport regime in the classical sense. However, it may so happen that there are systems where certain assumptions may need to be relaxed, since the ordering of various moments and their derivatives may become different due to unique physical conditions of that system. In this section, we consider two special cases of interest and outline the assumptions which are relaxed. However we do not solve the equations here, as the solution method remains the same as the previous sections, and merely adds some extra terms to the linearized system of Eqs ( 44) and (45). And so, we proceed to describe the two special cases, the first being the extension of the closure to account for the terms present in the linearized Burnett approximation, from which we derive the second, i.e. the extension of the closure to account for certain comparable magnitude terms in the drift-approximation. A. 21N -moment multi-temperature closure consistent to the linearized Burnett approximation

In certain systems, it may so happen that the spatial gradients of the heat-flux h α and the stress tensor π α are found in the first-order of the Knudsen number K n , which corresponds to the linearized Burnett approximation [START_REF] Zhdanov | Transport processes in multicomponent plasma[END_REF][START_REF] Chapman | The mathematical theory of nonuniform gases[END_REF] . In such a case, the parallel 21Nmoment system gets modified as follows. Firstly, the equation for the heat-flux becomes

- 5 2
kT α m α ∇ .π α + 5 2 k m α n α kT α ∇ T α = β R 11 αβ , ( 56
)
where the parallel divergence (∇ .) = b( b.)(∇.), and the factor of 5/2 since

σ α = m α 2 c 2 α c α c α - 1 3 δc 2 α f α dc α - 7 2γ α π α , (57)
and since we continue to neglect π α .∇ T α . The equation for the stress tensor similarly becomes

4 5 5 2 n α kT α ∂w αr ∂x s + ∂h αr ∂x s + 2n α kT α ∂u r ∂x s = β R 20 αβ ,
which can be represented in terms of the usual heat flux q α as 4 5 direct solution through linear algebra techniques, and furthermore the expressions for rank-1 and rank-2 moments are coupled. Normally, they can be solved through the method of successive solutions, i.e. first determine a zeroth approximation h α and π α by ignoring the gradient and divergence terms in the 21N -moment scheme, then substitute these solutions into Eqs. (56) and Eq. ( 58) to obtain a first approximation, substitute the first approximation to obtain a second and so on until the solutions sufficiently converge. In principle, this requires one to perform the matrix inversion a sufficient number of times, and can be both algebraically tedious and computationally expensive. However, one can further assume that the higher order derivatives of the heatflux and the stress-tensor are higher than order one in K n , in case of which the first approximation will uncouple the equations, giving the solution expression for the heat-flux as

∂q αr ∂x s + 2n α kT α ∂u r ∂x s = β R 20 αβ , (58)
- 5 2
kT α m α ∇ • (η α 0 ) + 5 2 k m α n α kT α ∇ T α = β R 11 αβ , (60)
where η α 0 represents the zeroth approximation for the partial viscosity η α which we leave inside the divergence since it depends on the temperatures of the species. The parallel divergence term can be further expanded as

∇ .(η α 0 ) = β ∂η α 0 ∂T β (∇ T β ). + η α 0 ∇ • .
(61) The solution for the stress tensor becomes

4 5 β ∇ λ βα 0 ∇ T β + α βα 0 w β + 2n α kT α = β R 20 αβ , (62)
where λ βα 0 is the zeroth solution for the partial parallel thermal conductivity and α βα 0 can be termed the zeroth solution to the partial parallel thermomotive coefficient respectively (i.e. the coefficient of the heatflux term depending on the diffusion velocities), which in turn depend on the temperatures of the species and can be expanded similarly as above in terms of gradients of temperature for the definition of the partial thermoelectric coefficients). These equations can now be solved individually rank-wise to obtain the first approximations to the partial viscosities η α 1 and partial thermal conductivities λ βα 1 (and partial thermoelectric coefficients similarly). The system of equations ( 60) and ( 62) may be considered the multi-temperature generalization of Eq. (5.2.21) of Ref. 7, and generally finds use in extended nonequilibrium thermodynamics, where the moment method is used for relating generalized thermodynamic forces to the thermodynamic fluxes [START_REF] Jou | Extended irreversible thermodynamics[END_REF][START_REF] Jou | Extended irreversible thermodynamics revisited (1988-98)[END_REF][START_REF] Zhdanov | The method of moments and the nonequilibrium thermodynamics of rarefied gases[END_REF][START_REF] Zhdanov | On a kinetic justification of the generalized nonequilibrium thermodynamics of multicomponent systems[END_REF] . Together with Eq. (59), they may be considered the extended 21Nmoment closure. On the basis of this, we now proceed to illustrate the closure in the drift approximation.

B. 21N -moment multi-temperature scheme in the drift-approximation

Generally, in the drift approximation of Hinton and Hazeltine [START_REF] Hinton | Theory of plasma transport in toroidal confinement systems[END_REF][START_REF] Helander | Collisional transport in magnetized plasmas[END_REF] , especially for tokamak regions with strong pressure gradients, it may so happen that the diamagnetic flow velocity u α,dia ∼ B×∇p α /n i Z i eB 2 +u α,dia is comparable to the main particle flow and the E×B flows. Accounting for this flow is important, since it often is the main contributor to the augmented friction force that leads to impurity accumulation in steady state with absence of turbulence [START_REF] Romanelli | Effects of density asymmetries on heavy impurity transport in a rotating tokamak plasma[END_REF][START_REF] Raghunathan | Heavy impurity confinement in hybrid operation scenario plasmas with a rotating 1/1 continuous mode[END_REF] . Furthermore, this diamagnetic flow is usually accompanied by the corresponding lowest order diamagnetic heat flux 39 q α,dia ∼ p α B × ∇ = (5/2)p α /n i Z i eB 2 + q α,dia , whose gradients may indeed be significantly larger than the gradients of particle fluid flow velocity and the E × B flow velocity. And indeed any further gradient of the gradient of the diamagnetic heat-flux would be an order smaller, and hence can be neglected. Thus, in such a situation, one needs to incorporate the diamagnetic heat-flux in Eqs. ( 56) and (58). In case of this, the expressions can be modified as follows to find

- 5 2 kT α m α ∇ •[η α 0 ( + α,dia )]+ 5 2 k m α n α kT α ∇ T α = β (R 11 αβ + M 110 αβ u α,dia ), (63)
and

4 5 {∇q α,dia } + 4 5 β ∇ λ βα 0 ∇ T β + α βα 0 w β + 2n α kT α ( + α,dia ) = β R 20 αβ , (64)
where α,dia = {∇u α,dia } . The equations for r α is similarly modified, but the equation for σ α remains the same, as follows

0 = β (R 12 αβ + M 120 αβ u α,dia ),
and 0 = β R 21 αβ .
Notice that there is no r α,dia in the rank-1 system of equations, since the value of q α,dia = (5/2)p α u α,dia is derived from setting h α,dia = 0, i.e. Grad's closure, and since r α,dia = (7/2γ α )h α,dia similarly. A partial solution for this system using transport coefficients obtained from the single-temperature Zhdanov closure is provided in Ref. 13, in which, the parallel divergence of η α 0 [ + dia, )] is neglected in Eq. (63), q α,dia is also neglected, and then a partial first solution is derived for the parallel divergence of π α , which is included in the momentum equation. It is also worth noting that including q α,dia in such a manner, either in the momentum equation or the closure or both, can lead to inconsistencies 40 such as summabilitity issues [START_REF] Poulsen | Turbulent transport in a multi-ion drift fluid model[END_REF] . That is, if a species is split into two continuous portions all else being equal, the two split portions may not evolve together as the whole unsplit species would. However, it becomes more and more consistent as the pressure profiles of various species approach each other. We can see that this form of the closure still remains quasi-linear (i.e.,at least linear in the derivatives), and it can be still solved through linear algebra techniques since the diamagnetic flow and its associated diamagnetic heat flux are fully determined (on knowing the pressure and density profiles and the magnetic field).

VII. SUMMARY OF CONCLUSIONS AND OUTLOOK

We first begin with the Sonine-Hermite polynomial moment method of Grad 42 and Zhdanov [START_REF] Zhdanov | Transport processes in multicomponent plasma[END_REF][START_REF] Alievskii | Transport equations for nonisothermal multispecies plasma[END_REF] , and using the methods described wherein, we re-derive and verify the most general moment-averaged balance equation presented in Appendix A of Ref. 7. We also explain certain subtleties of notation, leading to additional terms not present in Ref. 7. Using this general balance equation for the basis of linearization, we outline the assumptions behind the linearization of this system of balance equations, with detailed assumptions behind the order of terms, on the basis of which terms are retained and discarded. We note here generally that the plasmadynamical moments are considered of a zeroth order in Knudsen number K n , and the time derivative as well as the gradients of the plasmadynamical moments, and the higher-order moments are generally considered to be of first order in K n . The higher-order time derivatives and gradients of the plasmadynamical moments, and the first and higher-order time derivatives and the gradients of higher-order moments, are all considered to be of or smaller than the order (K n ) 2 .

On the basis of these approximations, we illustrate a general linearization scheme (following, but not exactly the same as, Secs. 4.6 and 6.3 of Ref. 7), which leads to a hierarchy of higher-order parallel moment equations such that the only equations with a non-vanishing LHS are the ones for the parallel heat-flux and the longitudinal stress tensor. The rest of the higher-order moment equations only contain the collisional RHS terms. This leads to the system of equations rendered solvable by linear algebra techniques. The solution expresses the higher-order parallel vectorial moments in terms of the flow velocities and temperature gradients, and express the higher-order longitudinal tensorial moments in terms of the rate-ofstrain tensor, which resemble the classical transport solutions in their usual form. In particular, the coefficient of the parallel temperature gradients in the parallel heatfluxes are the partial parallel thermal conductivities, and the coefficient of the longitudinal rate-of-strain tensor in the longitudinal stress tensor is the partial longitudinal viscosity. Furthermore, we illustrated the calculation of the generalized friction and thermal forces using the expressions for the higher-order moments found in this manner, which is what is key for implementation in code packages which simulate SOL/edge plasmas. Thus, having illustrated the general Zhdanov closure, we proceeded to study its convergence with respect to addition of a number of extra moments. To quantify this a little better for different impurities found in fusion, in addition to deuterium-tritium, we chose Carbon and Argon, mid-weight impurities at significant fraction of main ion density, and Tungsten, a heavy impurity at trace levels (Table I).

We choose seven rank-1 moments and eight rank-2 moments in order to test the convergence of the closure scheme. On calculating the total viscosity and the total thermal conductivities, the first observation we made was that the most significant addition to the viscosity and thermal conductivities seems to be indeed on going from the 13N -moment scheme to the 21N -moment scheme. The addition of further moments causes much smaller changes compared to the 21N -moment scheme. However, addition of more moments still leads to convergence of the scheme, particularly in representing the transport coefficients precisely to more decimal places. The total viscosity, in particular, seems to converge at a much faster rate than the total thermal conductivities. At this point, we chose an arbitrary two decimal places of precision to provide recommendations. We find again that the maximum differences for any transport coeffi-cient at any number of moments are found for the midweight impurities calculated with the single-temperature coefficients. This is especially pronounced for the thermal conductivities, on the basis of which we recommend using multi-temperature 21N -moment scheme for the mid-weight impurities to achieve the two decimal places precision. We furthermore observe that the differences between the single-temperature transport coefficients and the multi-temperature transport coefficients seem to still remain significant on addition of more moments, between 10 -80%, with the largest differences for the mid-weight impurities. Furthermore, we also calculate the augmentation of the friction and thermal forces caused by the additions of more moments. As with the convergence of transport coefficients, we find that the friction forces are more-or-less converged to two decimal places in the 21N -approximation, and the thermal forces to one decimal place or so. For the mid-weight impurities, particularly argon, the convergence for the thermal force was much slower. Furthermore, significant difference in the force values from the single-temperature and multi-temperature coefficients still persisted on addition of more moments. On the basis of this, we generally recommend using multi-temperature coefficients where feasible, in addition to the earlier recommendations.

Having found the 21N -moment scheme suitable for most cases, we proceeded to study the performance of the trace approximation against it. The trace approximation generally assumes the impurity particles to be present at trace values with respect to the background plasma density, and hence the trace approximation entails neglecting the impurity self-collisions and the back reaction of the impurities on the main plasma. On doing so, and on comparing the transport coefficients and forces, we found generally that the trace approximation tends to overestimate the transport coefficients and thermal forces by at least a factor of two for the mid-weight and low-weight impurities at significant densities. For the trace tungsten case, we found that the transport coefficients and the forces practically overlapped, justifying the trace limit. Generally, we also found that there persisted significant differences between the single-temperature and multitemperature cases as well (at least 40 %). On the basis of these, we can recommend using the multi-temperature trace approximation on the 21N -moment Zhdanov closure for heavy impurities at trace levels in order to speed up computations in fluid codes.

At last, we also described methods by which some of the assumptions may be relaxed, particularly two cases, the first in which the closure needs to respect the linearized Burnett approximation, and following which, the closure in the drift approximation. We did not compute the transport coefficients for these relaxed cases, as essentially the matricial method remains the same, but just adds more relatively constant terms.

In the scope of this article, we studied closures in the multi-species linear transport regime to obtain parallel transport coefficients, but we did not study the perpendicular ones. There are some analytical issues when a magnetic field is explicitly introduced in the multi-species system of equations in the linear transport regime, such as the increase in the degrees of freedom on which the transport parameters depend, as opposed to that of the ion-electron case where the degree of freedom essentially boils down to one variable β α = Ω α τ αα [START_REF] Ji | Closure and transport theory for highcollisionality electron-ion plasmas[END_REF][START_REF] Balescu | Transport processes in plasmas[END_REF] . The extension of the multi-species closure scheme to the fully magnetized case while addressing such issues, in order to study the perpendicular transport coefficients, is a subject of our future work. Though in the perpendicular direction, transport is often dominated by turbulence, they may still be useful to provide a more precise neoclassical description, which in recent literature has been found to account for impurity transport very well in absence of turbulence and MHD modes [START_REF] Angioni | Impurity transport in tokamak plasmas, theory, modelling and comparison with experiments[END_REF] . We also plan to consider the effects of parallel electric field on the closure and study the related thermoelectric effects. Furthermore, the scope of the study was committed to a fully analytical examination of the transport coefficients, and no concrete effects on numerical modelling were mentioned other than the suggested ranges of validity. The effects of the multi-species closure scheme with multitemperature coefficients on numerical SOL/edge fluid packages, Soledge3x-EIRENE in our case, is a subject we will explore.

Identities for the Sonine-Hermite polynomial G mn

The following relations for the Sonine-Hermite polynomials as given in Ref. 7, re-written in consistent terms, with derivation methods wherever needed 1. The derivative with respect to γ α becomes

∂G mn α ∂γ α = 1 2γ 2 α n(2m + 2n + 1)G m,n-1 α (A15) 2. Outer product with c α c αs G mn α = G m+1,n α + n γ α G m+1,n-1 α + 2 2m + 1 2m + 2n + 1 2γ α {G m-1,n α δ} s + {G m-
= n(2m + 2n + 1) 2γ 2 α × × G m+1,n-1 α + n -1 γ α G m+1,n-2 α + 2 2m + 1 2m + 2n -1 2γ α {G m-1,n-1 α δ} s + {G m-1,n α δ} s (A17)
This follows straightforwardly from the above two identities Eqs. (A15) and (A16).

∂c αr = nG m+2,n-1 α + n(n -1) γ α G m+2,n-2 α + 2n 2m + 3 δG mn α + n(2m + 2n + 1) (2m + 3)γ α δG m,n-1 α + {G mn αs δ} r - 2n 2m -1 δG mn αrs - n(2m + 2n + 1) (2m -1)γ α δG m,n-1 αrs + 2m + 2n + 1 4m 2 -1 × × 2{G m-2,n+1 α δδ} rs + 2m + 2n -1 γ α {G m-2,n α δδ} rs . (A19)
Notice the difference of the coefficient of δG mn αrs , and {G mn αs δ} r instead of G mn αs δ r from the expression in the appendix of Ref. 7. Notice the difference of a minus sign in 2m + 2n -1 in the coefficient of {G m-2,n α δδ}, note the specification of the rs th components in the δG m,n-1 αrs term. In the existing literature on symmetric and irreducible tensorial polynomials, there exist multiple notations, which often differ from one another by rational factors. This often makes the task of verifying expressions written in different notations difficult. Here we provide a short description of the three main notation systems encountered. The first form, notated without any special symbols, as used by Grad 42,[START_REF] Grad | Note on N-dimensional hermite polynomials[END_REF] and Zhdanov 7 write terms directly. For example, in this notation, ξδ is already symmetrized. For example

Double outer product with c

c αs c αr G mn α = G m+2,n α + 2n γ α G m+2,n-1 α + c 2 α 2m + 1 {δ r G mn α } s + c 2 α 2m + 1 {δ s G mn α } r + 1 4m 2 -1 (2m + 2n + 1)(2m + 2n -1) γ 2 α {G m-
{δ} ij = δ ij , {δ 2 } ijkl = δ ij δ kl + δ ik δ jl + δ il δ jk (B1) {ξδ} ijk = ξ i δ jk + ξ j δ ik + ξ k δ il (B2)
and so on. Some care has to be exercised with the notation, as in the example provided by Grad himself,

(ξ 2 + δ)δ = ξ 2 δ + 2δ 2 , (B3)
because δδ repeats symmetrized quantities in δ 2 twice. In such a notation, performing an outer multiplication of symmetrized irreducible rank-m tensor with ξ will need to be accompanied by a factor of 1/(m + 1) as such a multiplication would repeat the terms in the new rank-(m + 1) tensor (m + 1) times.

The second notation involves the use of parentheses (...) in the indices to indicate symmetrization. For example, we have the following

δ (ij) = δ ij , δ (ij δ kl) = 1 3 (δ ij δ kl + δ ik δ jl + δ il δ jk ) (B4) ξ (i δ jk) = 1 3 (ξ i δ jk + ξ j δ ik + ξ k δ il ) (B5)
For a general ξ (i1 ...ξ in δ j1j2 ...δ j2m-1j2m) , the relation between the two notations is given by

ξ n δ m = (n + 2m)! (2!) m m!n! ξ (i1 ...ξ in δ j1j2 ...δ j2m-1j2m) . (B6)
Proof of this relationship can be derived from simple combinatorics of indices. With this relation, one can verify that the Grad-Ikenberry polynomials [START_REF] Ikenberry | A System of Homogeneous Spherical Harmonics[END_REF][START_REF] Ikenberry | A system of homogeneous spherical harmonics[END_REF][START_REF] Ikenberry | Representation of Grad's Hermite polynomials as sums of products of sonine polynomials and solid spherical harmonics[END_REF] , the spherical tensor [START_REF] Weinert | Spherical tensor representation[END_REF][START_REF] Johnston | General Spherical Harmonic Tensors in the Boltzmann Equation[END_REF] and Grad's symmetric irreducible tensorial polynomials P (m) are one and the same, a fact often not clearly mentioned in previous literature.

The third notation encountered is the use of external curly brackets {...} on the terms being multiplied to indicate symmetrization. For example {ξ n δ m } indicates the symmetrization of the tensorial outer product of ξ n with δ m . This is essentially the same as the second notation ξ (i1 ...ξ in δ j1j2 ...δ j2m-1j2m) . It is sometimes found in quantum mechanics and continuum mechanics literature, but often just for notating simple symmetrization, e.g. of the sort when there are two tensors A and B of the rank-1, then the symmetric dyad {AB} = (AB + BA)/2. However, this is not so easy to understand for more complicated symmetrizations, where an index based notation is better.

In this article, we largely only use the first notation, and specify whenever it is used. The second notation and third notations are never used. However, a simple transformation can be performed in case required, using the abovementioned relationships between the different notations. However, we sparingly use curly brackets to denote symmetric, irreducible dyads, mentioning the definition when it is used. It is not to be confused with the third symmetrization notation.

Appendix C: Balance equations

One can derive balance equations for any given moment by substituting the appropriate values of (m, n) in Eq. (20). In this appendix, we derive a few specific cases for the purposes of verification of the LHS with known results, deriving general balance equations for scalar, vector and tensor moments, and finally for deriving the full balance equations for specific moments missing from previous literature.

13N -moment system of balance equations

The balance equation for the density n α b 00 α = ρ α is given by

dρ α dt + ρ α ∇ • u + ∂ ∂x r (ρ α w αr ) = β R 00 αβ . (C1)
Next, the balance equation for diffusion momentum

n α b 10 α = ρ α w α d dt (ρ α w α ) + ρ α w α ∇ • u + ∂ ∂x r (π αr + n α kT α ) -ρ α X α m α + Z α eE m α + u × ω α - du dt + ρ α w αs ∂u r ∂x s -w α × ω α = β R 10 αβ (C2)
Similarly, the balance equation for n α b 01 α = 0, which effectively results in a balance equation for the temper-ature T α is given by,

3 2 kn α dT α dt + 3 2 kn α w αr ∂T α ∂x r + ∂ ∂x r (h αr + n α kT α w αr ) - X αl m α + Z α eE αl m α + {u × ω α } l - du l dt ρ α w αl + π αrs ∂u r ∂x s + n α kT α ∇.u = β R 01 αβ , (C3)
from the fact that rst δ rs vanishes, and that rst π αrs ω αt vanishes because π αrs is symmetric. The curly brackets {AB} represent symmetric reduction of a dyad AB, e.g.

{A r B s } = 1 2 (A r B s + B r A s ) - 1 3 δ rs A l B l .
Notice that the evolution of temperature depends on the magnetic field only through the drift produced by the common flow. So far, the plasmadynamical equations take their usual form. The balance equation for the heat flux n α b 10 α = h α is calculated from

dh α dt + h α ∇ • u + ∂ ∂x r σ αr + 1 γ α π αr + 2 3 {n α b 02 α δ} r - X αl m α + Z α eE αl m α + {u × ω α } l - du l dt π αl -rst ω αt + ∂u r ∂x s n α b 30 αrs + 2 5 δh α + n α kT α δw α + {δh αs } r + 5n α k 2 w α dT α dt + 5k 2m α [π αr + n α kT α {δ} r ] ∂T α ∂x r = β R 11 αβ ,
where s, t, l are repeated indices summed over, and the free index r is used to show the remaining vectorial term from the double dot products. Now defining the third scalar moment n α b 02 α = θ α , the first rank-3 moment n α b 30 α = Ж α , and performing the expansions of δh α and δw α and their double dot products with ∂u r /∂x s , we get where π αsl rst ω αt can be thought of as the generalized tensorial cross-product (π αsl × ω α ) rs .

dh α dt + ∂ ∂x r σ αr + 1 γ α π αr + 2 3 ∇θ α -π αl X αl m α + Z α eE αl m α + {u × ω α } l - du l dt -h α × ω α +
The equations (C1)-(C5) may be considered as the most general moment equations derived from Grad's method in presence of electromagnetic and non-velocity dependent body forces. They agree with Eqs.(4.2.9)-(4.2.14) in Ref. 7. Notice that the most natural extension of this 13N -moment scheme would be to a 26N -moment scheme that incorporates the 1N additional moment(s) from θ α , 5N additional moments from σ α , and 7N additional moments from Ж α .

It is worth noting that this 13N -moment system of balance equations is not "regularized", meaning that it may not be applicable to cases involving propagation of shocks, doing which may introduce spurious entropy production [START_REF] Struchtrup | Regularized 13 moment equations for hard sphere molecules: Linear bulk equations[END_REF][START_REF] Rana | Thermodynamically admissible boundary conditions for the regularized 13 moment equations[END_REF][START_REF] Struchtrup | Macroscopic transport equations for rarefied gas flows[END_REF] .

General balance equations for moments with n ≥ 2 up to rank-2

The general balance equation for rank-0 scalar moments b 0n α , for n ≥ 2 is given by

d dt (n α b 0n α ) + n α b 0n α ∇ • u + ∂ ∂x r n α b 1n αr + n γ α n α b 1,n-1 αr -nn α b 1,n-1 αl X αl m α + Z α eE αl m α + {u × ω α } l - du l dt +n α ∂u r ∂x s -rst ω αt nb 2,n-1 αrs + n(n -1) γ α b 2,n-2 αrs + 2n 3 δb 0n α + n(2n + 1) 3γ α δb 0,n-1 α + n(2n + 1) 2
n α k m α b 0,n-1 α dT α dt + n(2n + 1) 2
n α k 2m α b 1,n-1 αr + n -1 γ α b 1,n-2 αr ∂T α ∂x r = β R 0n αβ . (C6)
The general balance equation for rank-1 vectorial moments b 1n , for n ≥ 2 is calculated as

d dt (n α b 1n α ) + n α b 1n α ∇ • u + ∂ ∂x r n α b 2,n αr + n γ α n α b 2,n-1 αr + 2 3 2n + 3 2γ α {n α b 0,n α δ} r + {n α b 0,n+1 α δ} r - X αl m α + Z α eE αl m α + {u × ω α } l - du l dt nn α b 2,n-1 αl + 2n + 3 3 {n α b 0,n α δ} l + n α ∂u r ∂x s -rst ω αt nb 3,n-1 αrs + n(n -1) γ α b 3,n-2 αrs + 2n 5 δb 1n α + n(2n + 3) 5γ α δb 1,n-1 α + {δb 1n αs } r + n α k m α n(2n + 3) 2 b 1,n-1 α dT α dt + n α k m α n(2n + 3) 2 b 2,n-1 αr + n -1 γ α b 2,n-2 αr + 2 3 2n + 1 2γ α {b 0,n-1 α δ} r + {b 0,n α δ} r ∂T α ∂x r = β R 1n αβ ,
which, on manipulating the terms as in the previous sub-section, can be written as In our previous article, we calculated the collision coefficients of partial bracket integrals in the form suggested by Chapman and Cowling 30 , Rat et al [START_REF] Rat | Transport properties in a two-temperature plasma: Theory and application[END_REF] . The calculation method in our previous article allowed for a free choice of a factor d αβ , such that on changing d αβ , one could find a new set of collision coefficients. For example, d αβ = 1 would represent a calculation similar to the original Chapman and Cowling calculation [START_REF] Bonnefoi | Contribution au calcul théorique des coefficients de transport d'un plasma d'azote par la méthode de Chapman-Enskog à l'approximation quatre de Sonine[END_REF][START_REF] Bonnefoi | Contribution à l'étude des méthodes de resolution de l'équation de Boltzmann dans un plasma à deux temperatures: Exemple le mélange Argon-Hydrogène[END_REF][START_REF] Rat | Transport properties in a two-temperature plasma: Theory and application[END_REF][START_REF] Chapman | The mathematical theory of nonuniform gases[END_REF] , d αβ = γ αβ /2 would represent the collision coefficients in a form similar to Zhdanov et al [START_REF] Zhdanov | Transport processes in multicomponent plasma[END_REF][START_REF] Alievskii | Transport equations for nonisothermal multispecies plasma[END_REF] , and we chose an additional d αβ = µ αβ /2kT as a reference. The coefficients were then compared to each other, in addition to a single-temperature coefficients, and a range of validity was provided for the the single-temperature coefficients. In this article, we restricted ourselves to only d αβ = γ αβ /2 for the multi-temperature coefficients as referring to only these as "multi-temperature" coefficients. We intend to demonstrate the reasoning for this here. There are some steps in the calculation that somewhat obscure a subtlety. The first is the definition of Π The exact definitions for the S, T , W's, X, g, a αβ , b αβ , X, ḡ can be found in Appendix A of our previous article [START_REF] Raghunathan | Generalized collisional fluid theory for multi-component, multi-temperature plasma using the linearized Boltzmann collision operator for scrape-off layer/edge applications[END_REF] . For the purposes of this section, it is sufficient to note that S and T are related to the exponential form of the Sonine polynomials, and that a αβ , b αβ , X, ḡ are expressions of S, T , W's, X, g. αβ as a convergent infinite series. However, the rate of convergence has not been addressed in the previous literature, as the series has not been expressed in this explicit form. This is the subtlety alluded to at the beginning of this note.

d dt (n α b 1n α )+ 2n + 5 5 n α b 1n α ∇•u+ ∂ ∂x r n α b 2,n αr + n γ α n α b 2,n-1 αr + 2 3 2n + 3 2γ α ∂ ∂x r (n α b 0,n α ) + ∂ ∂x r (n α b 0,n+1 α ) -nn α b 2,n-1 αl X αl m α + Z α eE αl m α + {u × ω α } l - du l dt - 2n + 3 3 n α b 0,n α X α m α + Z α eE m α + u × ω α - du dt -n α rst ω αt δ r b 1n αs +n α nb 3,n-1 αrs + n(n -1) γ α b 3
It is difficult to directly judge the rate of convergence of the series from this expression since it still depends on (S, T ) through b αβ . Therefore, it is instructive to try to see what this series aims to compute. We can define an H should reduce to H (m,p,q) αβ for any (p, q). Now, from the relations between the W's, X and ḡ, it is straightforward to write this as H (m,p,q) αβ ∝ exp -X 2 -b αβ ḡ2 S p m+1/2 (W 2 β )S q m+1/2 (W 2 α )

× (W β W α ) m P m ( Ŵ β . Ŵα )d X. (D11)
where X = X -(M α1 M α2 ) 1/2 (1 -θ αβ )ḡ, and b αβ = 1 -M α1 M α2 (1 -θ αβ ) 2 , and where P m is the Legendre polynomial. Now, if the Sonine and Legendre polynomials are series expanded, and then the W 's are expressed in terms of X 2 and ḡ( )2 , and integrate over 4π X2 d X, we will find that H (m,p,q) αβ evaluates to a finite series in ḡ2 and cos χ. The maximum power of ḡ2 in this series is p + q + m. Now, one can notice immediately that b αβ ḡ2 = [1 -M α1 M α2 (1 -θ αβ ) 2 ]ḡ 2 = (γ αβ /2)g 2 . Therefore, in the calculation of the bracket integrals in our previous article, for d αβ = γ αβ /2, the series for H (m) αβ converges perfectly on just retaining r ≤ p + q terms in the summation, with additional terms vanishing (because R m usually has a power of ḡ2m ).

It also implies that other choices of d αβ are not fully converged on retaining any finite number of terms, because in H which turns H (m,p,q) αβ into an infinite series. In our previous article, we had recommended a general limit of r ≤ p + q based on a comparison of terms expressed in terms of Chapman-Cowling integral and the series sum in H (m) αβ as was done in Refs. 25, 27, and 30. This explains why we do not observe singularities with the choice of d αβ = γ αβ /2, and that the singularities with other values of d αβ are a result of truncating the series too early. In Refs. 30 and 25, this leads to the correct result as b αβ = k αβ when the temperatures are equal, thus retaining only the r = 0 term. Thus, we take the opportunity to state that for the multi-temperature bracket integrals, only d αβ = γ αβ /2 will lead to a finite converged series for H (m) αβ with r ≤ p + q, and the rest will remain infinite sums as indicated in the expression. It may be possible to remedy the singularities by adding more terms to the series in the calculation of the bracket integral. Since however, a converged series is already available for d αβ = γ αβ /2, to do so would be a quixotic affair. Thus, the values of the collision coefficients generate with d αβ = γ αβ /2 may be considered the "closed forms" for the multi-temperature bracket integrals. Based on this, the generalization of the bracket

  + ∇.(n α c α G mn α ) -n α dG mn α dt + c α .∇G mn α + 1 m α F * α .∇ cα G mn α -c αs ∂G mn α ∂c αr ∂u r ∂x s = β R mn αβ .

  ) The coefficients A mnl αβ , B mnl αβ and C mnl αβ

7 FIG. 3 . 7 FIG. 4 .

 7374 FIG.3. Plot of augmentation percentage for friction forces for the multi-temperature (top) and the single-temperature (bottom) cases, on addition of moments. One can notice that the curves beyond the addition of the 13N -moment case are more-or-less overlapping. One can notice that the single-temperature coefficients do not suggest much variation in the augmentation over the temperature range.

FIG. 5 .

 5 FIG.5. Plot of friction and thermal force coefficients, i.e. the coefficients of the flow velocities and the temperature gradients in the RHS of the momentum equation, for the multi-temperature case, for M = 7, i.e. seven additional moments

λDFIG. 6 .

 6 FIG.6. Plot of total longitudinal viscosity η and total parallel thermal conductivities λα, for different species α,for multitemperature and single-temperature case, for the case of full collisions and trace approximation, for the 13N -moment and 21N -moment approximations, with respect to the temperature ratio. The curves the cool and warm colours represent the cases of full collisions and trace-approximation respectively.

with the reduced balance equations for n α b 12 α

 12 = r α and n α b 21 α = σ α remaining the same as earlier 0 = β R 12 αβ , and 0 = β R 21 αβ .

4 . 5 .

 45 Gradient with respect to c α Outer product of gradient w.r.t. c with c c αs ∂G mn α

--2n α {b 2n αsl rst ω αt } + n α k 2m α n(2n + 5 )b 2

 52 equation for rank-2 tensorial moments b 2n , for n ≥ 2, on similar manipulation of terms as for the previous rank-1 moments, is given bynn α b 3,n-1 αl X αl m α + Z α eE αl m α + {u × ω α } l -du l dt αβ . (C8)In general, one can also derive the moment-averaged balance equation for n < 2 from Eqs. (C6)-(C8), by simply setting the moments with negative indices to zero.Appendix D: The conundrum of d αβ : Why certain values don't work

  αβ = (1 -s) -5/2 (1 -t) -5/2 π -3 K αβ × {H (m) αβ (ḡ, χ) -H (m) αβ (ḡ, 0)}gσ αβ (g, χ)dΩdḡ, (D1)whereH (m) αβ (ḡ, χ) = exp -W 2 α -W 2 β -SW 2 β -T W 2 α × P (m) (W β ) : P (m) (W α )dX, (D2)by the absorption of the Sonine polynomials into the exponential. Then this is written asH (m) αβ (ḡ, χ) = exp -a αβ X2 -b αβ ḡ2 × P (m) (W β ) : P (m) (W α )d X. (D3)

0 ( 1 -

 01 The expression is then integrated over X, and what remains in the exponential is split and Taylor expanded as followsexp(-b αβ ḡ2 ) = exp(-k αβ ḡ2 ) exp((k αβ -b αβ )ḡ 2 ) = exp(-k αβ ḡ2 ) ∞ r=0 (k αβ -b αβ ) r r! ḡ2r , cos l χ)σ αβ (g,χ) sin χdχ, which are the d αβ values that are used in the current article. In the previous article, we had not given any recommendations for the choice of d αβ , only stating that the choice affected the form of the collision coefficients, and they made multi-temperature coefficients significantly different from the single-temperature ones when the temperatures are significantly different. This leads to a few peculiarities. For example, as long as the value of d αβ are positive and only depend on the masses and temperatures of the colliding species, one can have a multitude of different collision coefficients which in principle should be the same. However, once the collision potential is chosen, they provide values of collision coefficients which are different for different d αβ , seemingly making the solution multi-valued. Furthermore, in the current article, we find spurious singularities in the collision coefficients calculated from certain choices of d αβ . Now, with d αβ , the overall exponential term in H (m)αβ , on integrating over X can be written asH (m) αβ (ḡ, χ) ∝ exp -d αβ g 2 R m (S, T, ḡ, cos χ) × ∞ r=0 (d αβ g 2 -b αβ ḡ2 )where R m is some scalar function of ḡ and cos χ depending on the rank m of the tensor. It can be noted that, at this point this expresses H (m)

  (m,p,q) αβ of the following formH (m,p,q) αβ = exp -W 2 α -W 2 β S p m+1/2 (W 2 β )S q m+1/2 (W 2 α ) × P (m) (W β ) : P (m) (W α )dX. (D10)It is straightforward to notice that the expressions H (m) αβ

  (m,p,q) αβ with a different choice of d αβ , we would haveH (m,p,q) αβ ∝ exp -X 2 -k αβ ḡ2 × S p m+1/2 (W 2 β )S q m+1/2 (W 2 α )(W β W α ) m P m ( Ŵ β . Ŵα ) × ∞ r (b αβ -k αβ ) r r ! ḡ2r d X, (D12)

TABLE I .

 I Values of constants used for the different operational parameters. At 100eV, the maximum excitation state for higher-Z impurities is around +7, which the maximum charge states of Argon and Tungsten are limited to +7.

	α -β →	T-D		C-H			Ar-H	W-H
		nα	10 19		10 18			10 17	10 14
		Zα	+1		+6			+7	+7
		mα	3 amu		12 amu	40 amu	184 amu
		Tα	100 eV		100 eV	100 eV	100 eV
		n β	10 19		10 19			10 19	10 19
		Z β	+1		+1			+1	+1
		m β	2 amu		1 amu	1 amu	1 amu
	Total Viscosity η	0.02 0.04 0.06 0.08 T 0.6 0.00	0.8	1.0	1.2 Viscosities	1.4	1.6	1.8	2.0	1 2 3 4 5 6 7 8
				Temperature ratio TT/TD

β 50 -200 eV 50 -200 eV 50 -200 eV 50 -200 eV FIG. 1. Plot of longitudinal viscosity w.r.t. the temperature ratio for multi-temperature d αβ = γ αβ /2 case on addition of multiple moments.

Each number M in the legend represents, in principle, an additional 5M -moments. One can immediately notice that only the addition of σα as a moment contributes significantly, and every moment after that makes a much smaller contribution, so small in fact that the lines mostly seem to overlap beyond M = 2.

  Each number in the legend represents, in principle, an addition of 3M -moments. One can generally observe the same sort of convergence scheme as for the total longitudinal viscosity for the coefficients λD and λT . Notice however the relatively slower convergence of the electron counterpart λe.

				Thermal Conductivity λ e
		15
									1
									2
		10							3
	λe								4
									5
									6
									7
		0.6	0.8	1.0	1.2	1.4	1.6	1.8	2.0
				Temperature ratio TT/TD
				Thermal Conductivity λ T
		220
		200							1
		180							2
									3
	λT	160							4
		140							5
									6
		120
									7
		100
		0.6	0.8	1.0	1.2	1.4	1.6	1.8	2.0
				Temperature ratio TT/TD
				Thermal Conductivity λ D
		700
		600							1
		500							2
	λD	400							3 4
		300							5
		200							6
		100							7
		0
		0.6	0.8	1.0	1.2	1.4	1.6	1.8	2.0
				Temperature ratio TT/TD
	future work. At this point, we have been able to verify
	that similar patterns follow for the generalized Zhdanov
	closure as observed previously in closures of Braginskii's
	and Balescu's type.

[START_REF] Braginskii | Transport Processes in a Plasma[END_REF]

(since they are all weakly proportional to nkT τ and nk 2 T τ /m respectively, where τ is the mean collision time -see Ref. 15 for details). The scaling between λ e and λ α (α = D, H) is also of the order of (m e /m α ) 1/2 as per Braginskii's estimations as well. A full comparison of further convergence and comparison with Braginskii's closure is a part of the planned FIG.

2

. Plot of total parallel thermal conductivities λα, for different species α,for multi-temperature d αβ = γ αβ /2 case with the addition of more moments. The topmost plot is for λe, the middle one for λD and the bottom one for λT .

TABLE III .

 III .396 15.547 47.759 48.639 48.421 67.327 29.996 27.274 65.231 14.453 0.997 74.049 2 30.892 15.876 48.073 48.679 48.423 66.654 30.302 28.538 62.008 13.723 1.001 72.570 3 30.800 15.884 48.094 48.679 48.424 66.636 30.903 29.098 61.453 13.372 1.005 71.864 4 30.785 15.882 48.084 48.680 48.424 66.637 31.060 29.389 61.315 13.376 1.005 71.871 5 30.787 15.881 48.079 48.681 48.424 66.637 31.171 29.550 61.218 13.427 1.005 71.972 6 30.786 15.880 48.076 48.683 48.425 66.638 31.222 29.647 61.132 13.449 1.005 72.017 7 30.786 15.879 48.076 48.685 48.426 66.639 31.257 29.708 61.066 13.451 1.004 72.022 Differences in the augmentation of the friction force on addition of more moments. .425 68.190 5.283 0.819 52.278 15.449 15.318 40.846 17.836 15.233 41.451 3 44.665 7.122 68.051 5.738 0.808 53.029 16.560 15.761 46.712 13.880 13.899 48.659 4 44.350 7.100 68.382 5.769 0.808 53.080 17.037 16.193 47.433 14.104 13.448 49.461 5 44.465 7.099 68.248 5.768 0.806 53.082 17.315 16.563 47.283 13.792 13.403 48.856 6 44.443 7.099 68.273 5.768 0.804 53.084 17.588 16.810 47.289 13.869 13.424 48.524 7 44.442 7.099 68.275 5.766 0.804 53.081 17.808 16.969 47.407 13.852 13.441 48.471

	M	D-T			H-C			H-Ar			H-W
	I	II	III	I	II	III	I	II	III	I	II	III
	1 30M	D-T			H-C			H-Ar			H-W
	I	II	III	I	II	III	I	II	III	I	II	III
	2 44.953 7

TABLE IV .

 IV Differences in the maximum augmentation of the thermal force on addition of more moments.

  Civita tensor rst is antisymmetric to switching of indices. Without loss of generality, one may also write the term in the above equation n α kT α w αs ∂u r ∂x s + n α kT α w αt ∂u t ∂x r + n α kT α w α ∇.u kT α {w α δ} r + {h α δ} r

			= 2n α kT α w αs	∂u r ∂x s	+	5 3	n α kT α w α ∇.u,
	as it is found in some literature.
		And finally, the balance equation for the stress tensor
	n α b 20 α = π α is calculated as
		dπ α dt n α -+π α ∇•u+ ∂ ∂x r Ж αr + 2 5 5 2 X αl m α + Z α eE αl m α + {u × ω α } l -du l dt	ρ α {w α δ} l
	+	∂u r ∂x s -dπ α dt +π α ∇•u+∇.Ж α +	4 5	5 2	n α kT α		∂w αr ∂x s	+	∂h αr ∂x s
		-ρ α w αr	X αs m α	+	Z α eE αs m α	+ {u × ω α } s -	du s dt
		+ 2 π αrl	∂u s ∂x l	+ 2n α kT α		∂u r ∂x s		-2{π αsl rst ω αt }
												=	R 20 αβ , (C5)
													β
																Ж αrs	∂u r ∂x s
																+	7 5	h αs	∂u r ∂x s	+	2 5	h αt	∂u t ∂x r	+	7 5	h α ∇.u
																+ n α kT α w αs	∂u r ∂x s	+ n α kT α w αt	∂u t ∂x r	+ n α kT α w α ∇.u
																+	5 2	n α kw α	dT α dt	+	5 2	k m α	π αr	∂T α ∂x r	+	5 2	k m α	n α kT α	∂T α ∂x r	=

β R 11 αβ ,

(C4)

where rst Ж αrsl ω αt vanishes because Ж αrsl is symmetric to the switching of any two indices and the Levi-

rst ω αt [{δπ αs } r + n α kT α {δδ} rs ] = β R 20 αβ

On manipulating some terms, without loss of generality, and on realizing that rst {δδ} rs = 0, this takes the form

ACKNOWLEDGMENTS

The projects leading to this publication have received funding from Excellence Initiative of Aix-Marseille Université -A*MIDEX, a French Investissement d'Avenir Programme, project TOP & AMX-19-IET-013.

DATA AVAILIBILITY

The data in the article can be made available by the authors on reasonable request.

Appendix A: Polynomial Identities

We we collect together some properties for the polynomials provided by Zhdanov [START_REF] Zhdanov | Transport processes in multicomponent plasma[END_REF] and Weinert et al [START_REF] Weinert | Spherical tensor representation[END_REF] , rewritten in our notation.

Identities for Sonine polynomials S n m

The Sonine polynomials S n m (x) follow the recurrence relations

for any scalar x, The derivative with respect to x is given by

The gradient w.r.t. c α is given by

And finally, the derivative with respect to γ α is given by

2. Identities for the Irreducible tensorial monomial P (m)

For the irreducible monomial P (m) (ξ), the following relations, all in the symmetrization notation, all but the last from Ref. 7 1. The Rodrigues form of the irreducible tensorial monomial P (m) (ξ) is given by

where n!! indicates the double factorial of n. Sometimes, the symbol for the spherical tensor Y m (ξ) is used in place of P (m) (ξ), as they are both the same. (See Appendix B) for why the spherical tensor is equivalent to the irreducible monomial.)

2. Contraction with ξ r , i.e. 1-fold inner product with ξ, ξ • P (m) a1a2...ar...am ,

This is essentially the same as the contraction identity of the spherical tensor Y m ≡ P (m) in Ref. 45.

Gradient ∇P

where δ r refers to symmetrization without the participation of the r th index. E.g. ξδ r = ξ s δ rt +ξ t δ rs , and δ r δ s = δ ru δ st + δ rt δ su .

Additionally, one can also define derivative with respect to γ, if ξ = γ 1/2 c,

4. Outer product with ξ, ξP (m) ,

or

This identity can also be verified as in Refs. 45-47, through using the Rodrigues form of the spherical tensor directly.

Tensorial double gradient

rs .

(A12)

Note the difference from the identity provided in the appendix of Ref. 7.

6. Outer product of the gradient with ξ, ξ∇P (m) ,

For all the above identities, any r th or (rs) th component of a tensor mentioned as P refers to the additional ranks the tensor has, i.e. the rank-0 and rank-1 quantities do not possess an additional r th or (rs) th component, and rank-2 quantities do not possess an additional (rs) th component. In case such identities call for terms which have P (0 or 1) r or rs or P (2) rs , such terms can be set to zero. integral calculation to an arbitrary rank-m is going to comprise a part of our future work.