Equations and improved coefficients for paralleltransport in multicomponent collisional plasmas: method and application for tokamak modelling

New analytical expressions for parallel transport coefficients in multicomponent collisional plasmas are presented in this paper. They are improved versions of the expressions written in [V. M. Zhdanov. Transport Processes in Multicomponent Plasma, vol. 44. 10 2002.], based on Grad's 21N-moment method. Both explicit and approximate approaches for transport coefficients calculation are considered. Accurate application of this closure for the Braginskii transport equations is discussed. Viscosity dependence on the heat flux is taken into account. Improved expressions are implemented into the SOLPS-ITER code and tested for deuterium and neon ITER cases. Some typos found in [V. M. Zhdanov. Transport Processes in Multicomponent Plasma, vol. 44. 10 2002.] are corrected.


Introduction
Fusion toroidal devices with magnetic confinement (tokamaks and stellarators) operate with multispecies plasma. Plasma composition can contain mixtures of: main components: hydrogen, deuterium, tritium [2] and helium isotopes ; and impurities: helium, lithium, beryllium, carbon, nitrogen, neon, argon, etc. [3]. (In some operational regimes, helium might be the main component and then, in addition to the usual impurities, isotopes of hydrogen would be impurities.) Edge plasmas usually are in a collisional regime, by which we mean that macroscopic parameters change slowly on parallel (perpendicular) length scales of the mean free path (gyroradius) and time scales between collisions (gyromotion period). Under this assumption, a closure method such as the one proposed by Braginskii [4] can be applied to the moments of the distribution function. However, this approach is applied for a single ion species case and assumes only trace levels of impurities. For the multicomponent case, Grad's method [5] can be used. This method is based on the tensorial Hermite polynomials finite expansion approximation of the distribution function with a local Maxwellian distribution function as the zeroth-order approximation. It is developed in detail in [1]. Based on it, the 21N-moment approximation is studied in this work, though only parallel components are considered in this paper. Equations for electrons can be solved separately due to the large electron-ion mass difference, as it was considered in [1]. For ions the situation is more complicated, especially for mixtures with comparable masses. This is analyzed in detail below.
Some previous work on implementing the moment approach can be found in [6] (based on [7,8] and described in [9]) where an explicit implementation (based on matrix inversions) was implemented in SOLPS; [10] where the multispecies closure was implemented into 3D turbulence and transport codes based on [1]; and additional work in [11]. Some of this work assumed m light /m heavy 1 and so could not be applied to deuterium, tritium and helium mixture cases, while in others the closure for viscosity was not considered. In addition, there is a difference between definitions made by Braginskii [4] and in [1]. Thus, for the application of the closures discussed in [1] for the Braginskii equations, special corrections are developed in this paper (see Section 2). Some typos in [1] which directly affect the results are corrected in appendix C. In [6], a full closure (including viscosity) was implemented into the transport code, though the heat flux dependent part of viscosity was not considered and this higher order effect plays an important role in toroidal systems [12,13]. Such viscosity was added by Rozhansky et al [14,15] for correct calculation of the radial electric field in the H-mode pedestal in tokamaks, but only in the single fluid case and can be extended to multispecies cases.
In this paper this generalization of the equation system is presented. Section 2 is dedicated to an accurate application of Grad's closure discussed in [1] for the Braginskii system of equations used in transport codes, for instance SOLPS-ITER [16,17]. In addition, the heat flux dependent part of the viscosity is considered in section 2.
Since the calculation of transport coefficients is supposed to be performed during the numerical solution of algebraic equations, this might be a problem in a fluid code there this would be done at each time-step in each grid cell; therefore new improved analytical expressions are developed in section 3. In the case where masses of species are significantly different, such methods can be applied. This method is an improvement of the method discussed in [1]. Comparison for the improved expressions, the original expressions (8.4.7) in [1], and the solution of the explicit system of algebraic equations is provided. As it is shown there, even for the cases, where masses of components are not too different, analytical formulae provide surprisingly good agreement to the explicit numerical approach. Moreover, for some cases, the new formulae provide a much better match with the numerical results than the original analytical expressions (8.4.7) in [1]. Results can be found in appendix A.
Then, in section 4, a comparison between transport coefficients calculated using Grad's 21N-moment closure and transport coefficients used in the current SOLPS-ITER model is considered. Additionally, new expressions for thermal and friction forces were implemented into the SOLPS-ITER code and tested for neon transport in an ITER deuterium plasma and compared with the current SOLPS-ITER model [18,19].
Details can be found in [4]. In this equation system temperatures, heat fluxes and viscosities are defined according to Braginskii: where f αZ is the distribution function for ion or electron species αZ, k and l are component indices, and integrals of the collisional term C αZ due to Coulomb (elastic) collisions provide: Terms S n αZ , S m αZi , S E αZ describe, respectively, particle, momentum, and energy sources due to inelastic collisions between charged species (ionization, recombination, excitation) and all interactions with neutrals, and they are considered as external parameters in this model. The momentum flux due to the flow velocity is: The density and flow velocity are defined in appendix A.
In a multispecies plasma, the transport coefficients, which help to express q αZ through densities, velocities, temperatures, and their spatial derivatives, differ significantly for different plasma species, different species densities and different mass ratios of species. This means that no simplification may be made in the general case, and general methods for the solution of kinetic equations should be used.
Such a method is Grad's method of 21N-moment, described in [1]. Following [1], we define: where: One can note that temperature (10), heat flux (11) and viscosity (12) are defined in [1] regarding to the mass-averaged flow velocity (13), while (4)-(6) are defined with respect to the flow velocity of corresponding species. Definitions (10)- (12) are more suitable for this method, therefore they were used in [1]. Thus, the result of the closure discussed in [1], that is expressed in heat flux, viscosity, friction term and heat exchange term, can be applied for Braginskii system of equations (1)-(3) using corrections due to definition difference. Indeed: A similar correction should be applied to the heat exchange term: where according to [1]: where summation is performed over all species and time between collisions definition can be found in appendix A. One can recognize a heat source due to friction between different ions in the second term of (17).
The friction term R αZ defined according to (7) does not need correction provided that Coulomb collisions do not lead to particle sources and sinks and three possible definitions are equivalent: Therefore the friction term R αZ found using the approach discussed in [1] may be directly substituted into the system (1) -(3).
Note, according to Eq. (8.1.3) in [1], that momentum conservation in collisions is maintained: So is energy conservation according to (17): Details about the conservative property of collisions can be found, for instance, in [4].
Thus, assuming that the 21N-moment method is applied and transport coefficients are found, we need then to switch from h αZ and ← → π αZ to q (Br) αZ and ← → π (Br) αZ by using (15), (16) and taking into account the connection between temperatures (14).
In the present paper we consider only so-called parallel (with regard to the Bfield) transport coefficients. The classical transport across the magnetic field is usually of little consequence for magnetic fusion devices, since the anomalous transport is in most cases much larger. Then, as it is shown in [1], the transport coefficients appearing in the heat flux and friction force may be calculated independently of the viscosity coefficients. Therefore we may consider their corresponding calculations separately.
Finally, this approach can be applied both for electrons and ions. However, due to the small electron-ion mass ratio, transport for electrons can be considered separately, which significantly simplifies the approach. It is discussed in detail in [1]. Consequently, all the analysis in this paper will be devoted to ion transport coefficients.

Heat flux and friction term
Applying the Grad method with 21N-moment (see [1]), one can express h αZ through the velocities w αZ and temperature gradients ∇T αZ with the help of kinetic coefficients that can be found solving the algebraic system equations (8.4.2) in [1]. For the parallel (with regard to the B-field) component of h αZ , one gets (see Eq. (8.4.6) of [1])): where ∇ = (b · ∇), and summation is performed over all types of ions and: Collisional times τ (Zh) βα and average squared charge Z 2 α can be found in appendix A. Kinetic coefficientsc are the result of solving the system of algebraic equations (8.4.2) in [1] and application of corrections for each charge state (details can be found in appendix A).
Note that the first two terms on the l.h.s. of (22) represent the thermal conductivity. The third term is the velocity difference dependent part of the heat flux that, in the Braginskii approach [4], only appears in the electron heat flux.
Then, following the method in [1], the friction term can be written in the same form as Eq. (8.4.5) of [1]: where coefficientsc are the result of solving the system of algebraic equations (8.4.2) in [1] and application of corrections for each charge state as well. In (24), the first two terms represent the thermal force, and the third term is the interspecies friction force.

Viscosity
Now let us consider the system of equations for viscosity. Note that, to describe parallel transport correctly [12, and references therein], it is required to take into account the heat flux dependence for the viscosity. However, this effect was not considered in the viscosity equations (8.1.6), (8.1.6') in [1]. The necessity to account for this heat viscosity requires a modification of the approach described in [1], namely, using the general expression for the moment equation (A1.8) in [1] and adding the heat flux dependent terms to the left-hand side of (8.1.6), (8.1.6') in [1] (after summation over charge states): where R 20 α is the collisional right-hand side of Eq. (8.1.6) in [1] summed over charge states, that depends on π β (β represents each ion in the mixture), and where R 21 α is the collisional right-hand side of Eq. (8.1.6') in [1] summed over charge states, that also depends on π β . W-tensors and collisional right-hand sides can be found in appendix A. The solution of this system of equations (25)- (26) can be expressed as: where dimensionless transport coefficientsc are the result of solving the system of algebraic equations (25) and (26), and the corrections for each can be found analytically, as it is done for the heat flux (appendix A). The first term in (27) is the velocity dependent part of the viscosity that is discussed in [1]. The last terms in (27) represent additional effects due to the heat flux that is taken into account in this paper by adding heat flux dependent terms into the left-hand side both in (25) and (26).
Usually, it is not a trivial task in complex geometry to take a divergence ∇· ← → π αZ , intended to be included into (2), even in case (27). One way to do it supposes to take into account only parallel and drift components of the velocity and parallel and diamagnetic components of the heat flux in the viscosity calculation. Details can be found in appendix B.
Finally, to apply the closure discussed in [1] and in this paper for the Braginskii system of equations (1)-(3), corrections of temperatures, heat fluxes and viscosities should be done according to (14), (15) and (16). Now, let us recall results of transport in the Pfirsch-Schlüter regime [12,13]. For a simple plasma, using the moment approach, parallel viscosity transport coefficients (in the Pfirsch-Schlüter regime) were found in [13]. In the single ion species case, one can calculate viscosity transport coefficients using (27) as well. Thus, in the variables used in [12] they are: That is close to what was obtained in [12,13] and applied in [14,15]. Therefore, in the collisional case, it is expected to obtain a radial electrical field close to the radial electrical field in the Pfirsch-Schlüter regime [20]. It is important to mention that to get µ α2 as in (28), it is necessary to add the heat flux dependent term both into (25) and (26). As a result, using this approach, the result (28) is extended to arbitrary plasma mixtures.

Summary
In this section we have considered an explicit method based on solving an algebraic system of equations, which can be applied for the various plasma compositions, for example deuterium, tritium, helium, and other impurities. This case occurs in current devices [21] and will be standard operating procedure in future reactors [22,23,24]. Thus it should be implemented into codes like SOLPS-ITER [16,17] that are used for fusion reactor operation predictions. Section 4 is dedicated to this.
In the next section, an analytical approach is considered that can be applied for many cases which is, however, less accurate for species with close masses.
3 Improved analytical expressions

Heat flux
Let us consider the closure in the case with one light and several heavy ion species. The resulting 2n s × 2n s matrix (where n s is the number of different types of species (which have different atomic nucleus) in a mixture) which is intended to be inverted to solve system of equations (8.4.2) [1], can be split into blocks: where each block can be written: Here and further index 0 represents the light species, while indices 1..n s correspond to heavy species. Indices q = 1, 2 and p = 1, 2 define corresponding blocks. Note that in the codes, which are applied for light main ion and heavy impurities case modeling [16], main and impurity species are usually explicitly specified, and equations for them are written differently. The approach discussed in this paper allow us to write equations uniformly for all species.
First of all, for the diagonal elements, estimations can be made (i = 0): Taking into account, that β=0..ns , cross elements in the first row and column are smaller than a qp 00 by at least the mass ratio. Thus,ã qp 00 elements of the inverted matrix can be obtained independently from the contributions due to heavy species in the cross elements andã qp 00 ã qp 0i . Consider the solution for the temperature dependent part of the heat flux (heat conductivity): whereã 11 αβ is an element from theÃ 11 part of the inverted matrix. In the most common case, where the density of the light species is higher than the density of the heavy species, the impact from cross elements on the light species heat flux is even smaller due to multiplication with n β . On the other hand, for the heavy species heat flux all terms in the sum are comparable; however it is easy to show that the heavy species heat flux in such case does not contribute significantly to the global heat balance and therefore T 0 temperature for light species averaged over all charge states. Temperature of the light species T 0 plays major role in the thermal force between light and heavy species (see subsection 3.2), and thus, determines heavy impurity transport.
Thus, in the case with one light and several heavy species, the transport coefficients can be obtained by solving equations for each type of species independently, ignoring cross elements. This allows us to derive analytical expressions for the transport coefficients (see appendix A).
These expressions are an improved form of the Zhdanov analytical expressions (8.4.7) in [1]. Zhdanov suggests obtaining transport coefficients keeping only terms with (m 0 /m i ) 0 order, which directly exclude cross elements from consideration (due to (31)), while the approach discussed in this paper keeps higher order terms (m 0 /m i ) n , where n = 1..∞ in diagonal elements. The reason for this action follows from the fact that the diagonal elements play a major role. Analysis of the matrix in the test cases confirms this fact: even for the mixtures, where masses of species become close, cross elements were smaller than diagonal elements. Therefore, higher order terms in diagonal elements affect the answer more than for the cross terms.
To prove this assumption, let us compare the heat conductivity for deuterium in the presence of impurity ions found using this approach with the result using explicit matrix inversion (for solving algebraic system equations (8.4.2) in [1]) and the Zhdanov analytical expression (8.4.7) in [1]. Consider the case (deuterium, carbon and another impurity), where the temperature gradient for both impurities is the same as the temperature gradient for deuterium, and the case, where the temperature gradient for both impurities is (for some reason) 2 times higher than for deuterium, to explore what role the difference in ion temperatures plays. In test cases considered in this paper, the amount of carbon and another impurity is chosen according to the rules Z 2 C n C /n D = 0.5, Figure 1: Temperature dependent part of the heat flux transport coefficient for deuterium using the explicit solution of the system of equations (8.4.2) in [1], the original Zhdanov formula (8.4.7) in [1] and the improved formula (appendix A) for the D + C + another impurity cases with equal distribution between charge states. The additional impurity is varied along the horizontal axis.
The dash-dotted line shows the case where: The number in parentheses at the top of the figure gives normalised impurity density (= Z 2 i n i /n D where Z 2 i is an averaged square of the charge of impurity defined in (50)) and is kept constant.
I n I /n D = 1.0, where "I" corresponds to another impurity and Z 2 i is the averaged square of the impurity charge defined in (50). The heat flux associated with the deuterium heat conductivity is: where the transport coefficient c is shown in Figure 1 and the collision time is defined in appendix A.
In Figure 1 different results are plotted: with similar masses and densities for light impurities and small ratios m D /m imp and n imp /n D for heavy impurities.
It is clearly shown that (m 0 /m i ) 0 order accuracy is not sufficient to get good agreement with the explicit matrix inversion result, while the addition of higher order terms only in the diagonal elements provides solutions close to the numerical method except in the hydrogen-helium region, where cross elements in the matrix become important. In this region, a different temperature gradient for the impurity affects the result because both close masses and densities in the impurity term in (33) play a role. The same comparison can be made for the velocity dependent part of the heat flux: where the transport coefficient c βα is shown in Figure 2. The velocity dependent part of the deuterium heat flux in the presence of high mass impurities predicted by both original (8.4.7) in [1] and improved analytical expressions is close to the result obtained from the matrix inversion approach. However, for impurities lighter than beryllium, the improved formula gives a solution much closer to the numerical one ( Figure 2) than the Zhdanov formula. On the other hand, for the high mass species, the calculated heat flux is less accurate (Figure 2). In the original monograph [1], it is suggested to set to zero transport coefficients that represent heavy-light species interactions in the heat flux (velocity dependent part) for heavy species. Although that the improved formula (74) in appendix A suppresses the coefficient by a factor (µ αβ /m α ) 3/2 for the cases, where m α > m β , numerical calculation provides even smaller velocity dependent part of the heat flux ( Figure 2). In the extreme case, where a difference between masses of species is large, for instance ion-electron mixture in the simple plasma, the velocity dependent part is set equal to zero for heavy species [1,4]. Thus, the reduction of the velocity dependent part for heavy species can be implemented by setting (74) equal to zero if m α /m β < 1. To keep all expressions the same for both low and high mass species, smooth analytical approximation to the step function can be applied.
Finally, following [1], the heat flux for each charge state can be found (appendix A).

Friction term
Now consider the friction term. According to (8.1.3) [1]: According to (36), one can conclude that α R α = 0 for any velocities, heat fluxes and r-moments (5-order vector moments see Eq. are symmetric regarding to the replacement of αβ with βα).
Note that the impact on the friction term of species α from the heat flux and rmoment of species β is proportional to µ αβ /m β and (µ αβ /m β ) 2 correspondingly. Hence, leading terms in both R 0 and R i , where i = 0, are dependent on the heat flux and r-moment of the light species and, therefore, on T 0 . Thus, where the heat flux/r-moment of heavy species are calculated inaccurately, it does not affect the friction term for heavy species (and light as well) summed over all charge states.
The friction term for each charge state can be found using: Substituting results of the heat flux and r-moment into (36) and (37), one can obtain analytical transport coefficients for the thermal and friction forces (appendix A).
Now we compare these results with the explicit matrix inversion method and the coefficients previously implemented into SOLPS-ITER [18,19], which are based on the Zhdanov analytical expressions (8.4.7) in [1].
First of all, note that, in the trace-impurity case (Z 2 i n i /Z 2 0 n 0 1, m 0 < m i ), transport coefficients in [18,19] for the friction force yield c (1) imp ≈ 1 and, for the thermal force, c (2) imp ≈ 1.56. However, according to (36), the thermal force depends on the masses of the participants (besides the factor √ µ αβ in the G-matrices). Indeed, formulas (99) and (94) in appendix A, using the notation from [18,19], provide: (here the second term in (94) for species 0 is absent, which becomes important if m 0 ≈ m i ). Thus, the light trace-impurity thermal force in [18,19] is larger than our improved expressions.
Then, consider the thermal force for each type of species summed over all charge states in the non-trace-impurity case for a case (deuterium, carbon and another impurity) where the temperature gradient for both impurities is the same as the temperature gradient for deuterium, and the case, where the temperature gradient for both impurities is (for some reason) 2 times larger than for deuterium, to find where difference in ions temperatures play a role. The transport coefficient for the thermal force averaged over all charge states for carbon and another impurity (Figures 3, 4) is described by: For a deuterium plasma with carbon, the original Zhdanov analytical formula results in a ∼ 1% deviation from the numerically calculated result. Therefore for carbon and heavier impurities this formula can be applied, although, a slight increase of deviation for the heavier impurity due to the presence of carbon is seen ( Figure 3). Thus, the thermal force averaged over all charge states for nitrogen and neon calculated using the 3.0.6 SOLPS-ITER version and higher (and as a result in [25,26,27,28] and others) is the same as can be obtained by numerically solving the system of equations for heat flux and r-moment (5-order vector moment) (8.4.2) in [1].
Furthermore, in the region of heavy impurities, the impurity temperature gradient does not play a role ( Figure 3). Indeed, the shape of the impurity distribu- However, in the light impurity region the situation is different. For cases where ∇ T D = ∇ T i , the analytical Zhdanov expression results in an up to 15% deviation in a (D+C+Li) mixture and up to 90% deviation in a (D+C+He) mixture ( Figure 4). If the temperature of the impurity is different, the deviation is even larger. In contrast, the improved analytical expression for the thermal force in (D+C+Li) gives only a 1% deviation for lithium ( Figure 4) and 6% deviation for carbon (Figure 3) from the numerically calculated result. And for the (D+C+He) case, the improved formula provides up to 10% deviation for helium ( Figure 4) and up to 16% for carbon ( Figure 3). Note that in this region, the temperature gradient of impurities becomes important.
Notice that deviations for carbon are due to the presence of significant amount of helium (n He /n D = 0.4) and lithium (n Li /n D = 0.2), for smaller amounts of helium/lithium this case is closer to the pure deuterium/carbon. In the (n He /n D = 0.04) case, the deviation from the explicit matrix inversion result is up to 5% for carbon, however becoming worse for helium: up to 30% (up to 50% for the case, where ∇ T He = 2 × ∇ T D ).
For the cases: (deuterium, carbon and hydrogen); and (deuterium, carbon and  [1] provides even the wrong sign (always positive). However, the new formula (94) in appendix A works surprisingly well, which confirms the observation made in the tests that even for comparable masses, cross elements play a minor role. In mixtures (D+C+H/T). the thermal force has a deviation for hydrogen/deuterium/tritium of up to 30% (up to 50% for the case, where ∇ T T = 2 × ∇ T D ) and for carbon 15-20% deviation from numerically calculated result (figures 3, 4). For more accurate calculations in case of comparable masses, the explicit solution of the system of equations (8.4.2) in [1] is required. Moreover, for such cases the temperature gradient of impurities plays a role, and heat fluxes of impurities are important and an accurate calculation of these fluxes is also necessary. Now consider the thermal force for each charge state species obtained using (37), which affects space separation between different charge states of the same type of species. Note that this procedure does not require additional assumptions, thus can be made straightforwardly even for mixtures of species with close masses (appendix A). What was done for the thermal force in [18,19] represents what is done in the first term in (37), therefore the dependence on the difference between heat flux and r-moment (5-order vector moment, see Eq. (8.1.2) in [1]) averaged over charge states and for each charge state was not taken into account.
Here again consider the case where the temperature gradient for impurities is the same as the temperature gradient for deuterium, and the case, where temperature gradient for impurities (all charge states) is 2 times larger than for deuterium, to find, where the difference in ion temperatures plays a role. In Figure 5, the transport coefficient for the thermal force for the first charge state for carbon and another impurity described by: is plotted.
According to the thermal force for each charge state, expression (97) in appendix A, that is written taking into account the correction to the average thermal force, the first charge state is affected by the correction more than the others, because the impact from the average thermal force is reduced by a factor 1/Z 2 α . It is clearly seen ( Figure 5) that, for the first charge state of the impurity, the correction, represented by the third and fourth term in (37) (blue curves represent the model from [18,19]), plays a role. As a result, the temperature gradient for impurities becomes important ( Figure 5). However, for carbon and higher mass impurities, the thermal force that leads to space separation of each charge state, is an order of magnitude smaller than the thermal force that drives impurities, summed over all charge states ( Figure 3). Therefore, it is expected that for the heavy impurities, their temperature gradient affects slightly the the charge states space distribution, but not the global impurity transport.
Then, consider the friction force for each type of species summed over all charge states. Take a look at the case (deuterium, carbon and another impurity). In Figure 6 are plotted transport coefficients for the friction force, averaged over all charge states, for carbon and another impurity, described by: where definitions can be found in appendix A.
The transport coefficient for the friction force between deuterium and impurity predicted by the original Zhdanov analytical expression (8.4.2) [1], that is used in [18,19], is underestimated significantly for all realistic cases compared to the direct numerical solution of the system of equations (8.4.2) in [1] ( Figure 6): 37% for helium, 23 % for carbon, 17% for neon. This difference in the friction force can affect impurity transport [29]. On the other hand, our improved expression is much closer to the numerical solution (for a discussion of the corresponding deviations see below). This is the result of the inclusion of next order terms in the mass ratio m 0 /m i . Indeed, one can see (appendix A) that, under assumption m 0 /m i = 0, our improved expression turns into the original Zhdanov formula (8.4.2) in [1].
Moreover, as it was discussed in the subsection 3.1 (Figure 2), for the heavy species, the velocity dependent part of the heat flux is not calculated precisely by our improved expression. The same result can be obtained for the velocity dependent part of the r-moment (5-order vector moment, see Eq. (8.1.2) in Figure 6: Velocity dependent part of the friction term transport coefficient for another impurity/deuterium (solid), for carbon/deuterium (dashed) and for another impurity/carbon (dotted) using the explicit solution of the system of equations (8.4.2) in [1], the original formula from [18,19], and our improved formula (appendix A) for the D + C + another impurity case. (See the caption to Figure 1 for additional information.) Figure 7: Zoom in the dotted box of Figure 6 [1]). This inaccuracy for the heat flux and r-moment affects the friction force. Therefore, it is necessary to suppress them to zero to get a better match with the numerical result. One way to do this is by setting coefficients (74) and (85) equal to zero, if m α /m β < 1. The result, an even better match (indicated as "Analytical MDF (this work)" in figures 6 and 7, while "Analytical (this work)" corresponds to unmodified heat flux and r-moment) with the numerical result, with the deviation from the matrix inversion solution: for tritium 5% versus 8%, for helium 2% versus 8%, for carbon 0.03% versus 0.7%, for neon 0.6% versus 2% (for carbon in presence of neon 1% versus 5%). However, the mismatch for another impurity/carbon friction force becomes worse ( Figure 6). Note that the transport coefficient, equal to unity for impurity/impurity interaction (used in [18,19]), is close to the explicit approach result for most cases. Besides, a correction for the friction force between different charge states can be applied (appendix A), though for heavy species this transport coefficient is ≈0.8-1.0 and significantly differs from unity only for light species.

Viscosity
Finally, using a similar approach for the viscosity equations (25) and (26), one can obtain analytical expressions for both the velocity and heat flux dependent parts of the viscosity (appendix A). Now again consider the deuterium, carbon and another impurity case. The viscosity for species α summed over all charge states (for comparison with the analytical result from the matrix inversion method, we assume W h β /n β = W hα /n α for arbitrary β): In figures 8 and 9, c are plotted. Other definitions can be found in appendix A.
It is shown (Figure 8) that, for light species, the original Zhdanov expression overestimates the explicit method result significantly, while the new formula provides a closer result. Close answers are obtained for the heat flux dependent part of viscosity using both the improved formula and matrix inversion calculation ( Figure 9). For this part the analytical expression in [1] was not presented. For the case m D /m imp ≈ 1 and n imp /n D ≈ 1, the analytical method works quite well (figures 8 and 9). However, for the heavy impurities, the analytical expression does not match the calculated answer, and this affects the impurity momentum equation. On the other hand, the total momentum equation (summed over all species) is mainly affected by the light component viscosity (in case m D /m imp < 1 and n imp /n D < 1), and µ α1 , µ α2 as discussed in [12] can be obtained in the Pfirsch-Schlüter regime in the presence of impurities using expressions from appendix A.

Summary
New analytical expressions discussed in this section and presented in appendix A can be applied for the cases with one light and several heavy ions, when heat flux and viscosity of the heavy species does not affect the solution (for example n imp /n D < 1). For plasmas where the mass of the main component is  not significantly different from the mass of the impurity, improved expansions provide solutions closer to the numerical one than can be obtained by the original expression in [1]. Moreover, for the lighter than main ions impurities, the sign of the new thermal force is inverted, therefore our new expressions provide a qualitatively correct result. Furthermore, new expressions are written in general form without specifying which species are main ions or impurities, and which species are light or heavy, as it was for original analytical expressions in [1,18,19]. Besides, these formulae describe transport for each charge state of all the ions. Therefore, it allows us to consider helium plasmas with impurities, which was not possible in the previous SOLPS-ITER model, where main ions were assumed to be singly ionized. However, either in the case of accurate transport calculation for mixtures with close masses or in the case where high order moments of the heavy species distribution function (heat flux, r-moment (5-order vector moment, see Eq. (8.1.2) in [1]), viscosity, σ-moment (4-order tensor moment see Eq. (8.1.2) in [1])) are important, the explicit approach discussed in the previous section should be applied. Finally, it is important to mention that these results can be applied only for the collisional plasmas, since time derivatives and gradients (which are of the order of (λ/L) 2 , where λ is a mean free path and L is a space scale of the gradients) were neglected in high order moments equations (heat flux, r-moment, viscosity, σ-moment) [1]. The only exception, the heat flux gradients are added into the viscosity and σ-moment equations (25)- (26) for taking into account the Pfirsch-Schlüter regime effects discussed in [12,13,20]. In case of the less collisional plasmas next order terms have to be considered, as well. ions are considered to be equal for all ions [15], therefore (3) has to be summed over all ions. It should be noted that even in this case, T αZ can be different due to (14). Second, the part of the heat flux that depends on the velocity difference between ions with different masses (last term in (22)), is not presented in models which are based on the Braginskii equations (where such a term appears only in the electron heat flux) and has been neglected in multicomponent cases [16,30,31,32,33]; in future this contribution should be added. Moreover, the effective ion conductivity also changes in plasmas with non-trace-impurities, even under the assumption of identical temperature for all ions. In Figure 10, it is shown how the ion heat conductivity calculated using the explicit matrix inversion approach be different from what is predicted by the current SOLPS-ITER model (section B.4) [34].
A comparison of the friction and thermal force between the current SOLPS-ITER model [18,19], our improved expressions, and the explicit approach is made in detail in section 3.2. Therefore, it is not repeated here. The improved formulae for the friction and thermal forces have been implemented into the SOLPS-ITER code. Test results can be found in section 4.3.
In the current SOLPS-ITER model (sections C.3.1 and C.4.1) in [34], the heat is written so as to ensure the global conservative property in collisions [4]. Using (17) for electrons and ions, the distribution between electron and ion channels can be found more accurately. One can use (17) to find the correct heat distribution between different ions, as well.

Viscosity
Consider an orthogonal curvilinear coordinate system where x is a poloidal coordinate and y is a radial coordinate (43), as it is done in [14].
In this coordinate system (43) and under assumption that n αZ , T αZ , ϕ are flux surface functions inside the separatrix where parallel viscosity plays its larger role [12,13] in the parallel momentum balance, Eq. (136) in appendix B turns into: Note that expression (44) is a generalization of the corresponding terms in the momentum equation in [14,15] for the multicomponent case. The new result are compared with that currently applied for SOLPS-ITER in Pfirsch-Schlüter regime [15] in Figure 10. Despite the fact that the effect on the viscosity (velocity part) π u ion summed over all ions is not more than 8% due to the new approach, the effect on the impurity viscosity (velocity part) is significant. For the high mass high-Z impurity (carbon and higher), the viscosity (velocity part) transport coefficient is ≈1.5-2.5 times higher than currently implemented in SOLPS-ITER. The heat flux dependent part of the viscosity π h ion (whose importance was discussed above) can be 20% higher for this closure (assuming for comparison W h β /n β = W hion /n for either β) than for the current SOLPS-ITER approach ( Figure 10). The impurity viscosity (heat flux dependent part) can be up to 1.9 times higher.

Test cases using improved expressions
In Figure 7 one can see that, even for deuterium plasmas with neon non-traceimpurity transport coefficient, the friction force calculated by the original Zhdanov formula Eq. (8.4.2) in [1] is different compared to our improved treatment, while the thermal force remains the same (Figure 3), that should affect neon transport in the tokamak. The new friction and thermal force treatments were implemented into the SOLPS-ITER code and tested for a deuterium and neon ITER case. We present two cases, a reference case, run with the original model for the friction and thermal force terms, and a simulation where our new treatment is applied. These cases are cataloged in the ITER Integrated Modelling Analysis Suite (IMAS) database as #123077 and #123078, respectively. For the test, an intermediate case between cases 1b (#123014) and 2b (#123018) (with drifts) was chosen from those that are discussed in [25], with parameters: divertor neutral pressure p n = 7.5P a and relative Ne concentration, separatrixaveraged, c N e = 1.0%. To illustrate the changes in the impurity transport due to the new thermal and friction force treatments, the neon density is plotted in Figure 11.
One can see that, even for the deuterium and neon plasma, the mass difference is not high enough to assume m D /m N e ≈ 0. The application of the new formulas can observably affect the impurity transport for such cases. Also, after  Figure 11: 2D neon density distribution for the two identical ITER modeling cases using the SOLPS-ITER code, in which the original expressions from [18,19] are applied (#123077) (a) and our improved expressions are applied (#123078) (b) .
implementation of the new expressions, the relative separatrix-averaged Ne concentration dropped to c N e = 0.8%, while the neon seeding rate is kept constant (6.0 · 10 19 particles/s). In these tests, neon throughput and pumping speed below the Dome are kept constant. However, actual extraction of neon from the divertor by the pumping system, which depends on the distribution between neon flows in the divertor (part goes into the pumping system and part goes upstream), slightly increases, that was seen on the run diagnostics during the simulation. This leads to 10% drop of the amount of neon in the whole computational domain and in the core in particular. It is interesting to note that the amount of neon in the outer divertor region decreases by 20%, while the amount of neon in the inner divertor region does not change significantly. This impact on the neon distribution should be investigated.
It is also interesting to note that in the regions with low amounts of neon, the thermal force coefficients are slightly smaller, while the friction force remains the same according to (38), which is different compared to the non-trace-impurity regions, where transport coefficients follow behaviors described in figures 3 and 7. This example again shows that an accurate treatment of the main component distribution function shape affects the impurity transport both for the traceimpurity and non-trace-impurity cases. We demonstrate these runs to show that our improved expressions for the transport coefficients affect impurity transport in tokamaks like ITER. Further detailed analysis is required to understand this impurity behaviour and is supposed to be conducted in the following papers.
It is necessary to mention that, for this test, the expressions for thermal (97) and friction forces (106) in appendix A were slightly reduced. The thermal force for the impurity type "I" and charge state "Z" is assumed to be: The friction force for the impurity type "I" and charge state "Z" is assumed as: where "0" is connected to the main light single charge state ion and thermal and friction forces for the main ions are: On the one hand, such approach allows us to change only transport coefficients for this first test without rewriting equations significantly. On the other hand, for the deuterium and neon case these expressions provide close results to (97) and (106) in appendix A. However, full (97) and (106) will be implemented into SOLPS-ITER code later without specifying main ions and impurities explicitly.

Conclusions
This paper tackles the closure in the parallel direction of the system of fluid equations using Grad's method [5] for multicomponent collisional plasmas and extends the study made in [1]. The method allows one to obtain transport coefficients for the heat flux, friction and viscosity terms for each charge state of each species for the arbitrary plasma composition case, including deuterium + tritium + helium + other impurities mixtures, where densities and masses of species can be comparable. Therefore, this approach can be implemented into SOLPS-ITER [16,17] or other fluid codes to model this mixture for future reactors (ITER [22], DEMO [23], CFETR [24]) and complex compositions in existing machines.
In contrast to previous papers devoted to such an approach [10,6], two major improvements are made: corrections for the direct implementation of the closure [1] into the Braginskii equations; and the heat flux dependent part of viscosity. The importance of this part was discussed in [12,13]. In this paper the result in collisional regime is extended to arbitrary plasma mixtures.
New analytical expressions were developed (see appendix A) which provide a better match to the explicit matrix inversion method for the heat flux, friction and viscosity term and extend the applicability of the analytical approach to lower mass impurities. Moreover, this accurate treatment for multi charge state ions allows for the use of our analytical approach for cases where all ions have multiple charge states, for example helium plasmas with impurities, which was not possible with the original expressions.
These new friction and thermal force terms descriptions have been implemented into the SOLPS-ITER code and tested for an ITER deuterium plasma with nontrace-impurity neon and compared with the current SOLPS-ITER model [18,19]. Even for such a mixture, the new formulae show different impurity transport behavior confirming that the m D /m N e 1 assumption is not accurate enough. Thus, further studies of the impurity transport in ITER and other devices using this improved approach is required.
W-tensors in arbitrary Cartesian coordinate system are: where parallel velocity and drift contributions -diamagnetic and ExB drift velocities -should be taken into account in (58). Heat flux in (59) should contain the parallel contribution determined in this paper and the diamagnetic contribution (see appendix B).
For each charge state: where:
Averaged over all charge states: where: For each charge state: where:

A.2.3 Friction term
Friction terms can be split into thermal force and friction force components: Thermal force averaged over all charge states: where: which matches the corresponding (8.4.7) in [1].
For each charge state:
For each charge state: where:

A.3 S-coefficients
This section gives coefficients used for results in the previous section. They are used in [1] and [10]. These coefficients are obtained without additional assumptions, therefore they must match exactly the corresponding in [1] and [10].

B Appendix
This appendix considers the study of the divergence ∇ · ← → π (par) αZ , taking into account parallel and drift components of the velocity and parallel and diamagnetic components of the heat flux.
The parallel viscosity (27) in an arbitrary coordinate system is: where b = B/B and for u we substitute into W the sum of the parallel velocity, E × B and diamagnetic velocities: where summation is over all ions β with charge states ζ and For h α we substitute into W ha (using (51)), and for h αZ we substitute into W h αZ : where (see Eq.