© Q. Wargnier et al. 2020

1. Introduction

The lower atmosphere of the Sun is a complex and dynamic layer where the plasma is found in a wide range of different regimes – from weakly ionized and non-magnetized at the bottom of the photosphere to fully ionized and magnetized at the top of the transition region. In the Sun’s chromosphere, the pressure varies from a thousand pascals just above the photosphere to a few pascals in the high chromosphere (see Vernazza et al. 1981). Similarly, the magnitude of the magnetic field is large in active regions, around thousands of gauss in sunspots, whereas it is just a few gauss in quiet-Sun regions, as shown in Wiegelmann et al. (2014). It is still a great challenge nowadays to develop a unified model that can be used for both partially and fully ionized regimes under the large disparity of plasma parameters in the lower atmosphere.

The study of partially ionized plasmas in the presence of a magnetic field, such as in prominences and the lower atmosphere of the Sun, demands models that are beyond the ideal single-fluid magnetohydrodynamic (MHD) description. Phenomena necessary to fully understand the behavior of plasmas in the Sun chromosphere, such as Cowling’s resistivity, thermal conduction, heating due to ion-neutral friction, heat transfer due to collisions, charge exchange collisions, and ionization energy losses, are usually disregarded in ideal MHD models, or only described by means of ad-hoc terms.

Two types of fluid models for studying the lower Sun atmosphere are found in the literature. First, the single-fluid MHD description considers the plasma as a conducting fluid in the presence of a magnetic field. It has the main drawback of assuming thermal equilibrium conditions, where all the species are considered to be at the same temperature. This model is assumed to be valid at the photosphere, in a highly collisional framework, allowing us to study the formation of magnetic field concentrations at the solar surface in sunspots, magnetic pores, and the large-scale flow patterns associated with them (see Hartlep et al. 2012). It is also used for simulating the lower part of the atmosphere of the Sun, for example, incorporating subgrid-scale turbulence models for the transport of heat and electrical resistivity, as presented by Kitiashvili et al. (2015). The full MHD equations are solved in Martinez Sykora et al. (2015) accounting for non-grey radiative transfer and thermal conduction outside local thermodynamic equilibrium in order to study the effects of the partial ionization of the Sun chromosphere.

More recently, multi-fluid MHD models have been used to represent the non-equilibrium conditions of the chromosphere, based on continuity, momentum, and energy conservation equations for each species considered in the mixture (see Leake et al. 2013; Khomenko et al. 2014; Khomenko 2017; Alvarez Laguna et al. 2016, 2018; Ni et al. 2018). These models lead to very stiff systems that are difficult to solve numerically as they exhibit characteristic times that range from the convective and diffusive times of each fluid down to the collisional and chemical kinetics time scales, Alvarez Laguna et al. (2016). Leake et al. (2013) performed a multi-fluid simulation of magnetic reconnection for a weakly ionized reacting plasma, with a particular focus on the solar chromosphere, by considering collisional transport, chemical reactions between species, as well as radiative losses. Braginskii (1965) has derived rigorous expressions for the transport properties of fully-ionized plasmas starting from the Boltzmann equation. Khomenko et al. (2014) proposed a model for the description of a multi-component partially ionized solar plasma. Deriving rigorous transport properties for such multi-fluid model is complex and, so far, the theory has not yet been developed to the same level of accuracy as Braginskii’s.

In this paper, we propose a novel approach for studying the lower atmosphere of the Sun that is neither a single-fluid nor a multi-fluid MHD model, but a multi-component drift-diffusion model derived by Graille et al. (2009) from kinetic theory. This approach is able to capture most of the multi-fluid phenomena, that is different velocities between species, collisional exchange of mass, momentum and energy, chemical reactions, thermal non-equilibrium, and magnetized transport. This system of equations is less stiff as it solves for only one momentum equation, like in the single-fluid MHD approach. The electrons are assumed to move at non relativistic speeds in a collisional fluid regime. The transport properties are retrieved through a generalized Chapman-Enskog solution to the Boltzmann equation using a multiscale perturbation method. These developments lead to a model with an extended range of validity from partially to fully ionized plasmas, with or without the presence of a magnetic field. We couple this model to the Maxwell equations for the electromagnetic field and derive the corresponding generalized Ohm’s law for multicomponent plasmas. As in Braginskii’s theory, our model includes anisotropy in the transport properties of electrons that is created by the magnetic field. These properties are computed by solving for the integro-differential systems presented by Scoggins et al. (2016). We use a spectral Galerkin method based on the Laguerre-Sonine polynomial approximation previously studied in depth for various applications (see Devoto 1969; Ferziger & Kaper 1973; Woods 1995; Zhdanov 2002; Magin & Degrez 2004a; Bruno et al. 2011; Capitelli et al. 2013). The transport systems are implemented in the Mutation++ library that compiles state-of-the-art transport collision integral data for the different pairs of species in the mixture, see Scoggins & Magin (2014). In a simplified framework where the plasma is fully ionized, we compare the transport properties for the lower atmosphere of the Sun to the conventional expressions for magnetized plasmas attributed to Braginskii. For more general partially ionized conditions representative of the Sun’s lower atmosphere, we compute the muticomponent transport properties corresponding to the species diffusion velocities, heavy-particle and electron heat fluxes, and viscous stress tensor for a helium-hydrogen mixture in local thermodynamic equilibrium. Finally, the model is assessed for the 3D radiative MHD simulation of a pore in the highly turbulent upper layer of the solar convective zone. We compute the thermal conductivity, electrical conductivity, species diffusion coefficients, and components of the generalized Ohm’s law based on an existing fluid simulation. We draw conclusions about the importance of the contribution of its components, in particular, of the resistive and battery terms. The coupling of the model developed for transport with a flow solver is beyond the scope of this paper.

The structure of the paper is as follows. In Sect. 2, the non-dimensonal analysis used for the generalized Chapman-Enskog expansion is presented, together with the multi-component drift-diffusion model for two-temperature magnetized plasmas, the transport fluxes, and the generalized Ohm law. In Sect. 3, we describe the mixture considered, the conditions representative of the lower Sun atmosphere and the method used for computing the transport properties. In Sect. 4, we verify the model proposed on a fully-ionized case by comparing the results with those obtained by means of Braginskii’s theory. Finally, in Sect. 5, we discuss all the transport properties for a partially ionized case. Additionally, we compute the transport properties and the components of the generalized Ohm law for 3D radiative MHD simulations of a pore in the low Sun atmosphere.

2. Drift-diffusion model for multicomponent plasmas

In this section, we present the multi-component drift-diffusion model for two-temperature magnetized plasmas. It was derived from kinetic theory by Graille et al. (2009) as a generalized Chapman-Enskog solution to the Boltzmann equation, using a multi-scale perturbation method based on a dimensional analysis. Additionally, we compare this model to the multi-fluid description widely used for the Sun’s lower atmosphere.

2.1. Multi-scale analysis of the Boltzmann equation

We consider a multicomponent plasma composed of electrons, denoted here by the index 𝔢, and heavy particles (atoms and molecules, neutral or ionized), denoted by the subscript 𝔥. The species are assumed to be point particles, neglecting their internal energy. We combine the equations derived by Graille et al. (2009), Magin et al. (2009) for the fully magnetized case and the Maxwellian regime for reactive collisions.

In order to apply the Chapman-Enskog method, Graille et al. (2009) perform a multi-scale analysis on the non-dimensional Boltzmann equations for electrons and heavy particles. The order of magnitude of the different terms in the Boltzmann equation is studied through carefully chosen reference quantities. In the asymptotic fluid limit, the Knudsen number is assumed to be as the same order of magnitude as the square root of the mass ratio between electrons and heavy particles, defined as $\epsilon =\sqrt{m_{\mathfrak{e}}/m_{\mathfrak{h}}}$, where m_𝔢 and m_𝔥 are the mass of electrons and reference heavy particles respectively. This small parameter drives thermal nonequilibrium between the electron and heavy-particle baths. In the strongly magnetized regime, the Hall parameter is assumed to scale as ε⁰ = 1. The species distribution functions are expanded in the multiscale perturbation parameter ε following Enskog’s approach. As opposed to Braginskii, no assumption is made a priori on the zero-order distribution function. The asymptotic analysis of the Boltzmann equation is performed at successive orders of ε. The main results occurring at different time-scales are summarized in Table 1.

In time scales of order $t_{\mathfrak{e}}^{*}$, the electron population thermalize at the temperature T_𝔢. The electron distribution function is a Maxwell-Boltzmann distribution obtained by solving the electron Boltzmann equation at the order ε⁻². At order ε⁻¹ that corresponds to the time scale $t_{\mathfrak{e}}^{*}$, heavy particles thermalize at temperature T_𝔥. At the zeroth order ε⁰ that corresponds to the convective time scales, Euler equations for heavy particles and zero-order drift-diffusion equations for electrons are obtained. Finally, at order ε, corresponding to the diffusive time scale, we obtain Navier-Stokes equations for heavy particles and first-order drift-diffusion equations for electrons. It is important to mention that in Braginskii’s approach, the macroscopic equations are retrieved by taking moments of the Boltzmann equation by considering ad-hoc perturbations to the Maxwell-Boltzmann distributions. Only a correct scaling deduced from dimensional analysis can yield a sound multicomponent treatment.

2.2. Multi-component equations

We denote by symbol H the set of indices for the heavy particles of the mixture considered. First, the mass conservation equations for electrons and heavy particles are, as follows,

$$\begin{array}{c}\hfill {\partial }_t{\rho }_{\mathfrak{e}}+{\mathit{\partial }}_{\mathit{x}}·[{\rho }_{\mathfrak{e}}({\mathit{v}}_{\mathfrak{h}}+{\mathit{V}}_{\mathfrak{e}})]={\omega }_{\mathfrak{e}},\end{array}$$(1)

$$\begin{array}{c}\hfill {\partial }_t{\rho }_i+{\mathit{\partial }}_{\mathit{x}}·[{\rho }_i({\mathit{v}}_{\mathfrak{h}}+{\mathit{V}}_i)]={\omega }_i,i\in \mathrm{H}.\end{array}$$(2)

Here, ρ_𝔢 is the mass density of electrons, ρ_i is the mass density of heavy particle i ∈ H, v_𝔥 is the heavy-particle hydrodynamic velocity that has been chosen as the velocity reference frame, V_𝔢 is the electron diffusion velocity, V_i, i ∈ H is the diffusion velocity of heavy particle i in the reference frame, such that ∑_i ∈ Hρ_iV_i = 0. Quantities ω_𝔢 and ω_i, with i ∈ H, are respectively the chemical production rates of electrons and heavy particles.

Second, a momentum equation for all the particles within the plasma is found to be

$$\begin{array}{c}\hfill {\partial }_t({\rho }_{\mathfrak{h}}{\mathit{v}}_{\mathfrak{h}})+{\mathit{\partial }}_{\mathit{x}}·[{\rho }_{\mathfrak{h}}{\mathit{v}}_{\mathfrak{h}}\otimes {\mathit{v}}_{\mathfrak{h}}+p\mathbb{I}]=-{\mathit{\partial }}_{\mathit{x}}·{\mathbf{\Pi }}_{\mathfrak{h}}+\mathfrak{n}\mathfrak{q}\mathit{E}+\mathbf{I}\wedge \mathit{B},\end{array}$$(3)

where, p = p_𝔢 + p_𝔥 is the mixture pressure, which is the sum of the partial pressure of electrons p_𝔢 and heavy particles p_𝔥, Π_​𝔥 is the viscous stress tensor. The quantity ρ_𝔥 is the mass density of heavy particles, 𝔫 = 𝔫_𝔢 + 𝔫_𝔥 is the mixture number density (𝔫_𝔢 the number density of electrons and 𝔫_𝔥 the number density of heavy particles), 𝔫𝔮 is the mxiture charge defined by 𝔫𝔮 = 𝔫_𝔢𝔮_𝔢 + ∑_i ∈ H𝔫_i𝔮_i, where 𝔮_i is the charge of species i, E the electric field, and I the total current density defined as

$$\begin{array}{c}\hfill \mathbf{I}=\mathfrak{n}\mathfrak{q}{\mathit{v}}_{\mathfrak{h}}+{\mathit{J}}_{\mathfrak{e}}+{\mathit{J}}_{\mathfrak{h}}=\mathfrak{n}\mathfrak{q}{\mathit{v}}_{\mathfrak{h}}+n_{\mathfrak{e}}{\mathfrak{q}}_e{\mathit{V}}_{\mathfrak{e}}+{\displaystyle \sum _{i\in \mathrm{H}}}{\mathfrak{n}}_i{\mathfrak{q}}_i{\mathit{V}}_i,\end{array}$$(4)

where J_​𝔥 is the heavy-particle conduction current density, J_​𝔢 is the electron conduction current density, and B is the magnetic field. Note that electrons only participate in the momentum balance through the total pressure gradient and Lorentz force due to the large mass disparity between electrons and heavy particles.

Conservation of energy for free electrons and heavy particles is read as

$$\begin{array}{cc}& {\partial }_t({\rho }_{\mathfrak{e}}{e}_{\mathfrak{e}})+{\mathit{\partial }}_{\mathit{x}}·\left[{\rho }_{\mathfrak{e}}{e}_{\mathfrak{e}}{\mathit{v}}_{\mathfrak{h}}\right]+p_{\mathfrak{e}}{\mathit{\partial }}_{\mathit{x}}·{\mathit{v}}_{\mathfrak{h}}+{\mathit{\partial }}_{\mathit{x}}·{\mathit{q}}_{\mathfrak{e}}\hfill \\ \hfill & ={\mathit{J}}_{\mathfrak{e}}·{\mathit{E}}^{\prime }-\mathrm{\Delta }{E}_{\mathfrak{h}}^0-\mathrm{\Delta }{E}_{\mathfrak{h}}^1+{\mathrm{\Omega }}_{\mathfrak{e}},\hfill \end{array}$$(5)

$$\begin{array}{cc}& {\partial }_t({\rho }_{\mathfrak{h}}{e}_{\mathfrak{h}})+{\mathit{\partial }}_{\mathit{x}}·\left[{\rho }_{\mathfrak{h}}{e}_{\mathfrak{h}}{\mathit{v}}_{\mathfrak{h}}\right]+(p_{\mathfrak{h}}\mathbb{I}+{\mathbf{\Pi }}_{\mathfrak{h}}):{\mathit{\partial }}_{\mathit{x}}{\mathit{v}}_{\mathfrak{h}}+{\mathit{\partial }}_{\mathit{x}}·{\mathit{q}}_{\mathfrak{h}}\hfill \\ \hfill & ={\mathit{J}}_{\mathfrak{h}}·{\mathit{E}}^{\prime }+\mathrm{\Delta }{E}_{\mathfrak{h}}^0+\mathrm{\Delta }{E}_{\mathfrak{h}}^1+{\mathrm{\Omega }}_{\mathfrak{h}},\hfill \end{array}$$(6)

where ρ_𝔢e_𝔢 and ρ_𝔥e_𝔥 are the internal energies of electrons and heavy particles respectively, and q_𝔢 and q_𝔥, the electron and heavy-particle heat fluxes, respectively. Quantity E′ = E + v_𝔥∧B is the electric field in the heavy-particle reference frame, $\mathrm{\Delta }{E}_{\mathfrak{h}}^0$ and $\mathrm{\Delta }{E}_{\mathfrak{h}}^1$ are the energy relaxation terms at order ε⁰ and ε, respectively. Quantities J_​𝔥 ⋅ E′ and J_​𝔢 ⋅ E′ are the power that is developed by the heavy-particle and electron current densities, and Ω_𝔢 and Ω_𝔥 are the energy production rate for electrons and heavy particles.

Summing Eqs. (5) and (6), with an equation for kinetic energy, the total energy equation is obtained as

$$\begin{array}{c}\hfill {\partial }_t\mathcal{E}+{\mathit{\partial }}_{\mathit{x}}·\left[(\mathcal{E}+p){\mathit{v}}_{\mathfrak{h}}\right]+{\mathit{\partial }}_{\mathit{x}}·({\mathbf{\Pi }}_{\mathfrak{h}}·{\mathit{v}}_{\mathfrak{h}})+{\mathit{\partial }}_{\mathit{x}}·({\mathit{q}}_{\mathfrak{e}}+{\mathit{q}}_{\mathfrak{h}})=\mathbf{I}·\mathit{E},\end{array}$$(7)

where the total energy is $\mathcal{E}={\rho }_{\mathfrak{e}}{e}_{\mathfrak{e}}+{\rho }_{\mathfrak{h}}{e}_{\mathfrak{h}}+1/2{\rho }_{\mathfrak{h}}{{\mathit{v}}_{\mathfrak{h}}}^2$, and I ⋅ E the power developed by the electromagnetic field.

The system of Eqs. (1)–(7) is coupled to the set of Maxwell’s equations (8) :

$$\begin{array}{c}\hfill \begin{array}{cc}& {\mathit{\partial }}_{\mathit{x}}·\mathit{E}=\frac{\mathfrak{n}\mathfrak{q}}{{\epsilon }_0},\hfill \\ \hfill & {\mathit{\partial }}_{\mathit{x}}·\mathit{B}=0,\hfill \\ \hfill & {\partial }_t\mathit{B}=-{\mathit{\partial }}_{\mathit{x}}\wedge \mathit{E},\hfill \\ \hfill & {\mathit{\partial }}_{\mathit{x}}\wedge \mathit{B}={\mu }_0\mathbf{I}+{\mu }_0{\epsilon }_0{\partial }_t\mathit{E}\hfill \end{array}\end{array}$$(8)

where ε₀ is the vacuum permittivity and μ₀ the vacuum permeability.

The electron transport fluxes such as the electron diffusion velocity V_𝔢, the electron heat flux q_𝔢, the electron current density J_​𝔢 are composed of two terms: 1- a term at the convective time scale, at order ε⁰, corresponding to the Euler equations for heavy particles, and 2- a term which is a first order correction, at order ε, at the dissipative time scale corresponding to the Navier-Stokes equations for heavy particles. Similarly, the heavy-particle transport fluxes such as the heavy-particle diffusion velocity V_i, i ∈ H, heat flux q_𝔥, viscous stress tensor Π_​𝔥, and current density J_​𝔥 are defined at the dissipative time scale, at the order ε of the generalized Chapman-Enskog expansion.

The governing Eqs. (1)–(8) differ from the multi-fluid models used for partially ionized plasmas. Whereas multi-fluid models consider one hydrodynamic velocity distinct for each species, here only one common hydrodynamic velocity is used for the heavy species while each species diffuse in this reference frame. In addition, the structure of the governing equations is symmetrizable hyperbolic (Graille et al. 2009), which can be regarded as an important property for the numerical discretization of the system. Nevertheless, it is necessary to close the model by computing the transport properties. This computation is presented in Sect. 5 for a helium-hydrogen mixture. Additionally, by using the definition of the total current density I and the Maxwell equations Eq. (8), a generalized Ohm’s law for this particular model is derived in Sect. 2.5.

2.3. Transport fluxes for heavy particles

With the same formalism that is used by Graille et al. (2009), we introduce some extra notations in order to express the anisotropic transport properties in the presence of a magnetic field. First, a unit vector for the magnetic field ℬ = B/|B| is defined as well as the three direction matrices

$${\mathit{M}}^{\Vert }=\mathcal{B}\otimes \mathcal{B},{\mathit{M}}^{\perp }=\mathbb{I}-\mathcal{B}\otimes \mathcal{B},{\mathit{M}}^{\odot }=\left(\begin{array}{ccc}0& -{\mathcal{B}}_3& {\mathcal{B}}_2\\ {\mathcal{B}}_3& 0& -{\mathcal{B}}_1\\ -{\mathcal{B}}_2& {\mathcal{B}}_1& 0\end{array}\right),$$

such that we have for any vector x in three dimensions

$$\begin{array}{cc}& {\mathit{x}}^{\Vert }={\mathit{M}}^{\Vert }\mathit{x}=\mathit{x}·\mathcal{B}\mathcal{B},\hfill \\ \hfill & {\mathit{x}}^{\perp }={\mathit{M}}^{\perp }\mathit{x}=\mathit{x}-\mathit{x}·\mathcal{B}\mathcal{B},\hfill \\ \hfill & {\mathit{x}}^{\odot }={\mathit{M}}^{\odot }\mathit{x}=\mathcal{B}\wedge \mathit{x}.\hfill \end{array}$$

In the (x, ℬ) plane, the vector x^‖ is the component of x that is parallel to the magnetic field and x^⊥ is the perpendicular component. Therefore, we have x = x^‖ + x^⊥. The vector x^⊙ lies in the direction transverse to the (x, ℬ) plane. The three vectors x^‖, x^⊥, and x^⊙ are then mutually orthogonal. The anisotropic transport coefficients are expressed by means of the matrix notation

$$\begin{array}{c}\hfill \overline{\overline{\mu }}={\mu }^{\Vert }{\mathit{M}}^{\Vert }+{\mu }^{\perp }{\mathit{M}}^{\perp }+{\mu }^{\odot }{\mathit{M}}^{\odot }.\end{array}$$(9)

In the weakly or unmagnetized plasma regime the parallel and perpendicular component of the transport coefficients are identical, leading to μ^‖ = μ^⊥. In such framework, the transverse component of the transport coefficient (or in other words, perpendicular to the (x, ℬ) plane), is vanishing, leading to μ^⊙ = 0.

We recall that the heavy-particle transport fluxes presented in Sect. (2.2) are:

$$\begin{array}{c}\hfill {\mathbf{\Pi }}_{\mathfrak{h}},{\mathit{V}}_i,i\in \mathrm{H},\text{and}{\mathit{q}}_{\mathfrak{h}}.\end{array}$$

The viscous stress tensor Π_​𝔥 is defined as

$$\begin{array}{c}\hfill {\mathbf{\Pi }}_{\mathfrak{h}}=-{\eta }_{\mathfrak{h}}([{\mathit{\partial }}_{\mathit{x}}{\mathit{v}}_{\mathfrak{h}}+{({\mathit{\partial }}_{\mathit{x}}{\mathit{v}}_{\mathfrak{h}})}^{\top }]-\frac{2}{3}({\mathit{\partial }}_{\mathit{x}}·{\mathit{v}}_{\mathfrak{h}})\mathbb{I}),\end{array}$$(10)

where η_𝔥 is the viscosity of heavy particles. Then, the heavy-particle diffusion velocity V_i is given as

$$\begin{array}{c}\hfill {\mathit{V}}_i=-{\displaystyle \sum _{j\in \mathrm{H}}{D}_{\mathit{ij}}({\mathit{d}}_j+{\chi }_{\mathfrak{h}j}{\mathit{\partial }}_{\mathit{x}}ln{T}_{\mathfrak{h}}),i\in \mathrm{H},}\end{array}$$(11)

where D_ij is the multicomponent diffusion coefficient of heavy particles, d_j is the diffusion driving force of the particle j that is interacting with the heavy particle i, and χ_𝔥j is the heavy-particle thermal diffusion ratio of particle j. The diffusion driving force d_j is defined as

$$\begin{array}{c}\hfill {\mathit{d}}_j=\frac{1}{p_{\mathfrak{h}}}({\mathit{\partial }}_{\mathit{x}}{p}_j-{\mathfrak{n}}_j{\mathfrak{q}}_j{\mathit{E}}^{\prime }-{\mathfrak{n}}_j{\mathit{F}}_{j\mathfrak{e}}),j\in \mathrm{H},\end{array}$$(12)

which is composed of three forces: 1- the force due to the gradient of the partial pressure ∂_xp_j, 2- the Lorentz force, and 3- F_j𝔢 that is an average electron force acting on the heavy particle j. Quantity F_j𝔢 belongs to the category of diffusion driving forces and allows for a coupling between the heavy particles and electrons (Graille et al. 2009). This average force is defined as

$$\begin{array}{c}\hfill {\mathit{F}}_{j\mathfrak{e}}=-\frac{p_{\mathfrak{e}}}{{\mathfrak{n}}_j}{\overline{\overline{\alpha }}}_{\mathfrak{e}j}{\mathit{d}}_{\mathfrak{e}}-\frac{p_{\mathfrak{e}}}{{\mathfrak{n}}_j}{\overline{\overline{\chi }}}_{\mathfrak{e}j}{\mathit{\partial }}_{\mathit{x}}ln{T}_{\mathfrak{e}},j\in \mathrm{H}.\end{array}$$(13)

where ${\overline{\overline{\alpha }}}_{\mathfrak{e}j}$ and ${\overline{\overline{\chi }}}_{\mathfrak{e}j}$ are anisotropic transport coefficients. Finally, the heavy-particle heat flux reads

$$\begin{array}{c}\hfill {\mathit{q}}_{\mathfrak{h}}=-{\lambda }_{\mathfrak{h}}{\mathit{\partial }}_{\mathit{x}}{T}_{\mathfrak{h}}+p_{\mathfrak{h}}{\displaystyle \sum _{j\in \mathrm{H}}{\chi }_{\mathfrak{h}j}{\mathit{V}}_j+\sum _{j\in \mathrm{H}}{\rho }_j{h}_j{\mathit{V}}_j,}\end{array}$$(14)

where λ_𝔥 is the heavy-particle thermal conductivity and ρ_jh_j is the enthalpy of heavy particle j. The second term of Eq. (14) corresponds to thermal diffusion.

In the previous transport fluxes, some of the usual terms can be identified. The viscous stress tensor for heavy particles Eq. (10) is proportional to the strain tensor. Similarly, the first term of V_i in Eq. (11) is a generalized Fick’s law where the flux is proportional to the diffusion driving force d_j. Also, the first term of the heavy particle heat flux q_𝔥 in Eq. (14) is the usual Fourier law. In addition, V_i includes a term that is proportional to ∂_x​ln​T_𝔥, known as the Soret effect, described in Magin & Degrez (2004b), Giovangigli & Graille (2003, 2009).

In summary, the transport coefficients for heavy particles to be computed in the following sections are

$$\begin{array}{c}\hfill {\eta }_{\mathfrak{h}},D_{\mathit{ij}},{\chi }_{\mathfrak{h}j},{\lambda }_{\mathfrak{h}},i,j\in \mathrm{H}.\end{array}$$

In addition, the anisotropic transport coefficients associated to the coupling terms between electron and heavy particles are

$$\begin{array}{c}\hfill {\overline{\overline{\alpha }}}_{\mathfrak{e}j},{\overline{\overline{\chi }}}_{\mathfrak{e}j},j\in \mathrm{H}.\end{array}$$

2.4. Transport fluxes for electrons

The electron transport fluxes are

$$\begin{array}{c}\hfill {\mathit{V}}_{\mathfrak{e}}^{}\text{and}{\mathit{q}}_{\mathfrak{e}}.\end{array}$$

The electron diffusion velocity is defined as

$$\begin{array}{c}\hfill {\mathit{V}}_{\mathfrak{e}}=-{\overline{\overline{D}}}_{\mathfrak{e}}({\mathit{d}}_{\mathfrak{e}}+{\overline{\overline{\chi }}}_{\mathfrak{e}}{\mathit{\partial }}_{\mathit{x}}ln{T}_{\mathfrak{e}})+{\displaystyle \sum _{i\in \mathrm{H}}{\overline{\overline{\alpha }}}_{\mathfrak{e}i}{\mathit{V}}_i,}\end{array}$$(15)

where ${\overline{\overline{D}}}_{\mathfrak{e}}$ is the tensor for the diffusion coefficient of electrons, d_𝔢 is the electron diffusion driving force and ${\overline{\overline{\chi }}}_{\mathfrak{e}}$ is the electron thermal diffusion ratio. The electron diffusion velocity V_𝔢 is split into two terms : 1- the terms proportional to d_𝔢 and to ∂_x​ln​T_𝔢 are at order ε⁰, at the heavy particle convective timescale and 2- the terms proportional to V_i are at order ε at the heavy particle dissipative timescale. The electron diffusion driving force is defined as

$$\begin{array}{c}\hfill {\mathit{d}}_{\mathfrak{e}}=\frac{1}{p_{\mathfrak{e}}}({\mathit{\partial }}_{\mathit{x}}{p}_{\mathfrak{e}}-n_{\mathfrak{e}}{\mathfrak{q}}_e{\mathit{E}}^{\prime }),\end{array}$$(16)

which is composed of two forces: 1- the force due to the gradient of the partial pressure of electron ∂_xp_𝔢 and 2- the Lorentz force. The electron heat flux reads:

$$\begin{array}{cc}\hfill {\mathit{q}}_{\mathfrak{e}}& =-{\overline{\overline{\lambda }}}_{\mathfrak{e}}{\mathit{\partial }}_{\mathit{x}}{T}_{\mathfrak{e}}+(p_{\mathfrak{e}}{\overline{\overline{\chi }}}_{\mathfrak{e}}+{\rho }_{\mathfrak{e}}{h}_{\mathfrak{e}}){\mathit{V}}_{\mathfrak{e}}\hfill \\ \hfill & +p_{\mathfrak{e}}{\displaystyle \sum _{j\in \mathrm{H}}{\overline{\overline{\chi }}}_{\mathfrak{e}j}{\mathit{V}}_j+{\rho }_{\mathfrak{e}}{h}_{\mathfrak{e}}\sum _{j\in \mathrm{H}}{\overline{\overline{\alpha }}}_{\mathfrak{e}j}{\mathit{V}}_{\mathfrak{j}},}\hfill \end{array}$$(17)

where ${\overline{\overline{\lambda }}}_{\mathfrak{e}}$ is the electron thermal conductivity tensor and ρ_𝔢h_𝔢 is the enthalpy of electrons. Quantity q_𝔢 is split into two terms: 1- the terms proportional to ∂_xT_𝔢 and V_𝔢 are at the heavy-particle convective timescale, and 2- the terms proportional to V_i that are at the heavy-particle dissipative timescale.

As in the heavy species transport properties, some usual terms can be identified that are Fick’s and Fourier’s laws. Additionally, terms that are coupled to the heavy-particle diffusion are present at the first order of the generalized Chapman-Enskog expansion. In summary, the anisotropic transport coefficients associated to the transport fluxes for electrons are

$$\begin{array}{c}\hfill {\overline{\overline{D}}}_{\mathfrak{e}},{\overline{\overline{\chi }}}_{\mathfrak{e}},{\overline{\overline{\lambda }}}_{\mathfrak{e}}.\end{array}$$

In this section, a list of the transport fluxes and the corresponding transport coefficients has been presented. The method used for computing the latter will be presented in the next Sect. 3.

2.5. Generalized Ohm’s law

In the following, we derive a general expression for Ohm’s law in the previous set of governing Eqs. (1)–(8). In order to do so, we rewrite the expression of the electric current by grouping the terms in each of the driving forces. By doing this, we obtain a general algebraic expression for the electric field E as a function of the transport coefficients and the corresponding driving forces. Note that, in this derivation, the displacement current has been neglected.

With the electron diffusion velocity from Eq. (16) and the heavy-particle diffusion velocity from Eq. (11) we find the total current I as follows

$$\begin{array}{cc}\hfill \mathbf{I}& =\mathfrak{n}\mathfrak{q}{\mathit{v}}_{\mathfrak{h}}+\frac{{\left(n_{\mathfrak{e}}{\mathfrak{q}}_e\right)}^2}{p_{\mathfrak{e}}}{\overline{\overline{M}}}_{{\mathit{E}}^{\prime }}{\mathit{E}}^{\prime }-n_{\mathfrak{e}}{\mathfrak{q}}_e[{\overline{\overline{M}}}_{p_{\mathfrak{e}}}\frac{{\mathit{\partial }}_{\mathit{x}}{p}_{\mathfrak{e}}}{p_{\mathfrak{e}}}\hfill \\ \hfill & +{\displaystyle \sum _{j\in \mathrm{H}}{\overline{\overline{M}}}_{p_j}\frac{{\mathit{\partial }}_{\mathit{x}}{p}_j}{p_{\mathfrak{h}}}+{\overline{\overline{M}}}_{T_{\mathfrak{e}}}{\mathit{\partial }}_{\mathit{x}}ln{T}_{\mathfrak{e}}+{\overline{\overline{M}}}_{T_{\mathfrak{h}}}{\mathit{\partial }}_{\mathit{x}}ln{T}_{\mathfrak{h}}]}\hfill \end{array}$$(18)

where the multicomponent electromagnetic matrices $\overline{\overline{M}}$ are defined as:

$$\begin{array}{cc}\hfill {\overline{\overline{M}}}_{{\mathit{E}}^{\prime }}& =\frac{p_{\mathfrak{e}}}{p_{\mathfrak{h}}}{\displaystyle [\sum _{i\in \mathrm{H}}{\overline{\overline{\xi }}}_{\mathfrak{e}i}\left(\sum _{j\in \mathrm{H}}{D}_{\mathit{ij}}{\overline{\overline{\xi }}}_{\mathfrak{e}j}\right)+{\overline{\overline{D}}}_{\mathfrak{e}}],}\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill {\overline{\overline{M}}}_{p_{\mathfrak{e}}}& =\frac{p_{\mathfrak{e}}}{p_{\mathfrak{h}}}{\displaystyle [\sum _{i\in \mathrm{H}}{\overline{\overline{\xi }}}_{\mathfrak{e}i}\left(\sum _{j\in \mathrm{H}}{D}_{\mathit{ij}}{\overline{\overline{\alpha }}}_{\mathfrak{e}j}\right)+{\overline{\overline{D}}}_{\mathfrak{e}}],}\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill {\overline{\overline{M}}}_{p_j}& ={\displaystyle \sum _{i\in \mathrm{H}}{\overline{\overline{\xi }}}_{\mathfrak{e}i}{D}_{\mathit{ij}},j\in \mathrm{H},}\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill {\overline{\overline{M}}}_{T_{\mathfrak{e}}}& =\frac{p_{\mathfrak{e}}}{p_{\mathfrak{h}}}{\displaystyle [\sum _{i\in \mathrm{H}}{\overline{\overline{\xi }}}_{\mathfrak{e}i}\left(\sum _{j\in \mathrm{H}}{D}_{\mathit{ij}}{\overline{\overline{\chi }}}_{\mathfrak{e}j}\right)+{\overline{\overline{D}}}_{\mathfrak{e}}{\overline{\overline{\chi }}}_{\mathfrak{e}}],}\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill {\overline{\overline{M}}}_{T_{\mathfrak{h}}}& ={\displaystyle \left[\sum _{i\in \mathrm{H}}{\overline{\overline{\xi }}}_{\mathfrak{e}i}\left(\sum _{j\in \mathrm{H}}{D}_{\mathit{ij}}{\chi }_{\mathfrak{h}j}\right)\right],}\hfill \end{array}$$(23)

and the tensor ${\overline{\overline{\xi }}}_{\mathfrak{e}i}$ is defined as

$$\begin{array}{c}\hfill {\overline{\overline{\xi }}}_{\mathfrak{e}i}=\frac{{\mathfrak{n}}_i{\mathfrak{q}}_i}{{\mathfrak{n}}_{\mathfrak{e}}{\mathfrak{q}}_{\mathfrak{e}}}\mathbb{I}+{\overline{\overline{\alpha }}}_{\mathfrak{e}i},i\in \mathrm{H}.\end{array}$$(24)

After some algebra, using Eq. (18) and neglecting the displacement current, the general expression of the electric field is obtained as

$$\begin{array}{cc}\hfill {\mathit{E}}^{\prime }& ={\overline{\overline{M}}}_{{\mathit{E}}^{\prime }}^{-1}[\frac{p_{\mathfrak{e}}}{{(n_{\mathfrak{e}}{\mathfrak{q}}_e)}^2}({\mathit{J}}_{\mathfrak{e}}+{\mathit{J}}_{\mathfrak{h}})+\frac{p_{\mathfrak{e}}}{{\mathfrak{n}}_{\mathfrak{e}}{\mathfrak{q}}_{\mathfrak{e}}}({\overline{\overline{M}}}_{p_{\mathfrak{e}}}\frac{{\mathit{\partial }}_{\mathit{x}}{p}_{\mathfrak{e}}}{p_{\mathfrak{e}}}+{\displaystyle \sum _{j\in \mathrm{H}}}{\overline{\overline{M}}}_{p_j}\frac{{\mathit{\partial }}_{\mathit{x}}{p}_j}{p_{\mathfrak{h}}}\hfill \\ \hfill & +{\overline{\overline{M}}}_{T_{\mathfrak{e}}}{\mathit{\partial }}_{\mathit{x}}ln{T}_{\mathfrak{e}}+{\overline{\overline{M}}}_{T_{\mathfrak{h}}}{\mathit{\partial }}_{\mathit{x}}ln{T}_{\mathfrak{h}})].\hfill \end{array}$$(25)

The expression of the multicomponent electromagnetic matrices ${\overline{\overline{M}}}_{{\mathit{E}}^{\prime }},{\overline{\overline{M}}}_{p_{\mathfrak{e}}},{\overline{\overline{M}}}_{T_{\mathfrak{e}}},{\overline{\overline{M}}}_{T_{\mathfrak{h}}},{\overline{\overline{M}}}_{p_j}$ can be subdivided into two terms: (1) a term which depends on the coupled heavy-particle and electron transport properties, such as ${\overline{\overline{\alpha }}}_{\mathfrak{e}j},{\overline{\overline{\chi }}}_{\mathfrak{e}j},D_{\mathit{ij}},{\chi }_{\mathfrak{h}j}$, which scales at the dissipative timescale for the heavy particles at order ε, and (2) a term which depends only on the electron transport properties ${\overline{\overline{D}}}_{\mathfrak{e}},{\overline{\overline{\chi }}}_{\mathfrak{e}}$, which scales as the convective timescale for the heavy particles at order ε⁰.

Some usual terms can be identified in the general expression of the electric field Eq. (25). The first term of Eq. (25) is the resistive term, where the expression of the resistivity tensor is defined as

$$\begin{array}{c}\hfill {\overline{\overline{\eta }}}_{\mathfrak{e}}=\frac{p_{\mathfrak{e}}}{{(n_{\mathfrak{e}}{\mathfrak{q}}_e)}^2}{\overline{\overline{M}}}_{{\mathit{E}}^{\prime }}^{-1}.\end{array}$$(26)

The second term and third term of Eq. (25), can be identified as a general expression of the battery term for a multicomponent plasma due to the pressure gradients of electrons and heavy particles. The fourth and last term of Eq. (25) are additional terms due to the presence of Soret and Dufour terms in the equations of the diffusion velocities Eqs. (11) and (16).

In Appendix B, a simplified fully ionized plasma case has been considered which leads to a simplified expression for the electric field. In this case, the multicomponent electromagnetic matrices can be simplified, and the usual expression of the electric field and magnetic induction equation are retrieved (see in Appendix B).

3. Methodology

A helium-hydrogen mixture, composed of 92% Hydrogen and 8% of Helium in mole fractions, which is typical of the Sun’s atmosphere, Asplund et al. (2009), is considered. The set of species considered in this mixture is denoted as

$$\begin{array}{c}\hfill S_1=\{\mathrm{He},{\mathrm{He}}^{+},\mathrm{H},{\mathrm{H}}_2,{\mathrm{He}}^{++},{\mathrm{H}}^{+},e-\}.\end{array}$$(27)

The heavy species such as carbon, oxygen or metals are not taken into account. We assume that they do not impact the transport properties as they are trace elements, i.e., the mole fractions are very small.

We study the transport coefficients for the previous mixture within a range of temperature, pressure, and magnetic field that is largely representative of the lower Sun atmosphere, as shown by Vernazza et al. (1981), Carlsson & Stein (1995): the temperature varies from 1000 K to 30 000 K, the pressure from 1 Pa to 10⁴ Pa, and the magnetic field from a few Gauss to thousands of Gauss, as shown by Wiegelmann et al. (2014). In the following, the plasma beta parameter is defined as β_p = 2μ₀p/|B|², where p is the total pressure of the plasma in Pascal and |B| is the magnetic field in Tesla.

Consequently, for a range of temperature between 1000 K and 30 000 K, two cases have been considered. The case A, where the total pressure is p = 10⁴ Pa, and β_p = 10. The case B, where the total pressure is p = 1 Pa, and β_p = 0.1. Finally, the case A is chosen as a thermally pressure dominated case whereas, the case B is a magnetically pressure plasma case. In this framework, in the case A, the transport coefficients are shown to be isotropic and in the case B the latter are shown to be anisotropic. The latter example may be representative for conditions that we can find 1- in a cool sunspot, where the average temperature is generally around 4000 K and the magnitude of the magnetic field is about 1000 G, as shown in Fröhlich & Lean (2004), Solanki (2003), and 2- in the quiet sun photosphere, in the lower and upper part of the chromosphere, where the temperature is varying from 5000 K to 10 000 K, and the pressure from 10⁴ Pa to 0.1 Pa, as shown in Vernazza et al. (1981), Russell (1929), Asplund et al. (2009) (in orders of magnitude).

Based on the chosen conditions, we compute the thermochemical equilibrium composition. The mole fraction and the ionization degree of the helium-hydrogen mixture S₁ for the case A and case B are shown in Figs. 1–3. These results are obtained with a method that is based on the minimization of the Gibbs free energy with suitable mass balance constraints, developed by Scoggins & Magin (2015) in thermal equilibrium. The compositions that are shown in Figs. 1 and 2 will be used to study the transport properties in the following sections.

Fig. 1. Mole fraction of the helium-hydrogen mixture S1, for P = 104 Pa (case A), as a function of temperature.

Fig. 2. Mole fraction of the helium-hydrogen mixture S1, for P = 1 Pa (case A), as a function of temperature.

Fig. 3. Ionization degree the helium-hydrogen mixture of S1, for case A (full line) and case B (dashed line) as function of the temperature.

Similarly as Scoggins et al. (2016), Magin & Degrez (2004a), the calculation of the transport coefficients is based on the solution of integro-differential equations. In order to solve these equations, a spectral Galerkin method is applied. This method expands the coefficients in a series of orthogonal Laguerre-Sonine polynomials that are truncated at a given order of approximation. The calculation is thus reduced to a linear algebraic system of equations. As a result, the transport coefficients can be obtained by the resolution of linear systems with matrices that are known functions of the macroscopic parameters that are the field variables and the collision integrals between particles. The solution of these systems allows for the transport coefficients to be written as linear combinations of the collision integrals, which take into account the interaction potential for a collision between two particles. These linear combinations are derived by extending the definition and the calculation of bracket integrals introduced by Ferziger & Kaper (1973) or in Woods (1995), Zhdanov (2002), Balescu (1988) to the thermal nonequilibrium case, studied in depth by Kolesnikov (2003). According to Magin & Degrez (2004a), Kolesnikov (2003), Tirsky (1993), the transport coefficients involving collisions between heavy particles and electrons converge for expansions in second order non vanishing Laguerre-Sonine polynomials and higher. In this work, we use the third order Laguerre Sonine polynomials approximation in order to compute the transport properties. This method has been widely used in the literature, for example Bruno et al. (2011), Woods (1995), Magin & Degrez (2004b), Zhdanov (2002), Balescu (1988), Capitelli et al. (2013).

The explicit relations for ${\overline{\overline{\alpha }}}_{\mathfrak{e}j},{\overline{\overline{\chi }}}_{\mathfrak{e}j},$ and ${\overline{\overline{D}}}_{\mathfrak{e}},{\overline{\overline{\chi }}}_{\mathfrak{e}},{\overline{\overline{\lambda }}}_{\mathfrak{e}}$ in terms of the solutions to the transport systems can be found in Scoggins et al. (2016) and in Appendix C. The heavy particle transport systems for η_𝔥, D_ij, χ_𝔥j, λ_𝔥, are found in Magin & Degrez (2004b) with the difference that the mole fractions are given in terms of heavy species only, excluding electrons. These are also presented in Appendix C. The transport collision integrals for the helium-hydrogen mixture are taken from Bruno et al. (2010).

4. Verification of the method in a fully ionized plasma case S₂

In order to verify the presented method, we perform a comparison with Braginskii’s transport properties. In the method of Braginskii (1965), the computation of the transport properties as well as the derivation of the governing equations are valid only for fully ionized plasmas. As it can be seen in Fig. 3, the helium-hydrogen mixture S₁ can be considered to be fully-ionized, mainly composed of S₂ = {H⁺, e⁻}, when the temperature is higher than 15 000 K. The comparison will be thus performed in conditions where the mixture is S₂ in a range of temperatures from T = 15 000 K to T = 30 000 K for the case A and case B. To illustrate the comparison, we focus on the properties ${\lambda }_{\mathfrak{e}}^{\Vert }$, ${\lambda }_{\mathfrak{e}}^{\perp }$, η_𝔥 and λ_𝔥, although the rest of them show similar behaviour.

On the one hand, in Braginskii (1965), the derivation of the governing equations can be summarized in three main steps: 1- A fully ionized ion-electron plasma is considered in a constant magnetic field, 2- The Landau collision operators are used, simplified by the Lorentz process, and 3- an ad-hoc Chapman-Enskog method is used based on the square root of the mass ratio between electron and ions, Balescu (1988). On the other hand, in Graille et al. (2009), a general multicomponent plasma that can be partially or fully ionized is considered in a constant magnetic field, the Chapman and Cowling collision operators highly studied in Ferziger & Kaper (1973), Woods (1995) are used and the Chapman-Enskog expansion is performed after a dimensional analysis of the Boltzmann equation. Finally, the two methods lead to distinct governing equations.

Although the governing equations between the two models are different, the integro-differential systems for computing the transport properties are similar or even identical in the case of a fully ionized plasma. In both models, the anisotropic electron transport properties have the same integro-differential systems. However, only the systems related to the parallel component of the heavy particle transport properties are identical to those from the model derived by Graille et al. (2009). Consequently, only the parallel component of the heavy particle transport properties can be compared with those from the model of Graille et al. (2009). This is due to the fact that both models are based on the Chapman-Enskog expansion. However, the differences result from the scale analysis from the Boltzmann equation that is carried out by Graille et al. (2009) before applying the expansion.

In both models, the transport coefficients are expanded in a series of orthogonal Laguerre-Sonine polynomials. The latter are written as linear combinations of collision integrals that are simplified by potential interactions, based on the usual Coulomb interaction screened by the Debye-length. This approximation assumes collisions with large impact parameters and small scattering angles. However, in Braginskii (1965), the series are truncated at the second-order approximation whereas a third-order approximation has been performed in the case of the model presented. The expression of the transport coefficients depends on the mean collision times $\overline{{\tau }_{\mathfrak{e}}}$ and $\overline{{\tau }_{\mathfrak{h}}}$ defined as

$$\begin{array}{c}\hfill \overline{{\tau }_{\mathfrak{e}}}=\frac{3{m_{\mathfrak{e}}}^2{{\epsilon }_0}^2}{{\mathfrak{n}}_{\mathfrak{h}}{\mathfrak{q}}_{\mathfrak{e}}^{\mathfrak{4}}log(\mathrm{\Lambda })}{\left(\frac{2\pi k_B{T}_{\mathfrak{e}}}{{\mathfrak{q}}_{\mathfrak{e}}^{\mathfrak{2}}{\mathfrak{n}}_{\mathfrak{e}}}\right)}^{\frac{3}{2}},\overline{{\tau }_{\mathfrak{h}}}=\sqrt{\frac{2m_{\mathfrak{h}}}{m_{\mathfrak{e}}}}{\left(\frac{T_{\mathfrak{h}}}{T_{\mathfrak{e}}}\right)}^{\frac{3}{2}}{Z}^{-2}\overline{{\tau }_{\mathfrak{e}}},\end{array}$$(28)

where log(Λ) is the Coulomb logarithm defined by Spitzer (1963), k_B, Boltzmann’s constant, and Z, the charge number. The mean collision times as defined in Eq. (28), can be seen as a first-order Chapman Cowling approximation of the collision time for electron-ion, and, ion-ion collisions, as shown by Woods (1995). Correction terms depending on Z are used for the computation of the transport coefficients. This method leads to simplified expressions of the transport coefficients that depend only on the mean collision times and the charge number of the fully ionized plasma considered, as presented by Balescu (1988), Woods (1995).

In Braginskii (1965) (see Eq. (4.37)), the parallel and perpendicular components of the electron thermal conductivity tensor are defined as

$$\begin{array}{c}\hfill {\lambda }_{\mathfrak{e}}^{\Vert }\underset{\mathbf{Br}}{:=}\frac{{\mathfrak{n}}_{\mathfrak{e}}{\mathrm{k}}_{\mathrm{B}}^{\mathrm{2}}{\mathit{T}}_{\mathfrak{e}}}{m_{\mathfrak{e}}}\overline{{\tau }_{\mathfrak{e}}}\left[3.16\right],\end{array}$$(29)

$$\begin{array}{c}\hfill {\lambda }_{\mathfrak{e}}^{\perp }\underset{\mathbf{Br}}{:=}\frac{{\mathfrak{n}}_{\mathfrak{e}}{\mathrm{k}}_{\mathrm{B}}^{\mathrm{2}}{\mathit{T}}_{\mathfrak{e}}}{m_{\mathfrak{e}}}\overline{{\tau }_{\mathfrak{e}}}\left[\frac{4.664x^2+11.92}{x^4+14.79x^2+3.77}\right],\end{array}$$(30)

where Br denotes the computation of the transport coefficient as derived by Braginskii (1965). $x={\omega }_{\mathfrak{e}}\overline{{\tau }_{\mathfrak{e}}}$ and ω_𝔢 = 𝔮_𝔢B/m_𝔢 and the values in brackets correspond to Braginskii’s coefficients for a charge number Z = 1.

Figures 4 and 5 show the parallel and perpendicular component of the electron thermal conductivity tensor ${\overline{\overline{\lambda }}}_{\mathfrak{e}}$, as function of the temperature, for the case A and the case B, for the fully ionized plasma S₂. Here, we compare the expressions from Braginskii Eq. (29) and Eq. (30) with those that are given by Scoggins et al. (2016) that are based on a third-order Laguerre-Sonine polynomials approximation. Strong similarities are obtained in all the considered cases. In Braginskii (1965), the components of the electron thermal conductivity tensor are underestimated leading to differences that are less than 20%. These differences are increasing at high temperatures. Similar results have been obtained for all the other electron transport properties.

Fig. 4. Parallel component of the electron thermal conductivity tensor ${\lambda }_{\mathfrak{e}}^{\Vert }$ for a fully ionized plasma S2, as function of temperature. Dashed lines and full lines correspond to the transport coefficient from the model of Braginskii (1965), and from Graille et al. (2009) respectively. Bold lines correspond to the case A, the other lines correspond to the case B.

Fig. 5. Perpendicular component of the electron thermal conductivity tensor ${\lambda }_{\mathfrak{e}}^{\perp }$ for a fully ionized plasma S2, as function of temperature. Dashed lines and full lines correspond to the transport coefficient from the model of Braginskii (1965), and from Graille et al. (2009) respectively. Bold lines correspond to the case A, the other lines correspond to the case B.

Similarly, the parallel component of the heavy thermal conductivity and of the heavy particle viscosity of the model of Braginskii (1965, see Eq. (4.44)), have been compared with the expression from Scoggins et al. (2016) and Magin & Degrez (2004a). In Braginskii (1965), the heavy thermal conductivity and heavy viscosity are defined as

$$\begin{array}{c}\hfill {\lambda }_{\mathfrak{h}}^{\Vert }\underset{\mathbf{Br}}{:=}{n}_{\mathfrak{h}}{\text{k}}_{\text{B}}^2{T}_{\mathfrak{h}}\overline{{\tau }_{\mathfrak{h}}}\left[3.91\right],\end{array}$$(31)

$$\begin{array}{c}\hfill {\eta }_{\mathfrak{h}}^{\Vert }\underset{\mathbf{Br}}{:=}{n}_{\mathfrak{h}}{\text{k}}_B{T}_{\mathfrak{h}}\overline{{\tau }_{\mathfrak{h}}}\left[0.96\right].\end{array}$$(32)

Figure 6 and 7 show the heavy thermal conductivity λ_𝔥 and heavy viscosity η_𝔥 respectively, in the same conditions as in Figs. 4 and 5. As before, strong similarities have been obtained in all the considered cases for the chosen conditions, which leads to differences that are smaller than 20%. In addition, it can be shown that the heavy transport properties from Braginskii (1965) are isotropic at the chosen conditions.

Fig. 6. Heavy thermal conductivity λ𝔥 for a fully ionized plasma S2, as function of temperature. Dashed lines and full lines correspond to the transport coefficient from the model of Braginskii (1965), and from Graille et al. (2009) respectively. Bold lines correspond to the case A, the other lines correspond to the case B.

Fig. 7. Heavy particle viscosity η𝔥 for a fully ionized plasma S2, as function of temperature. Dashed lines and full lines correspond to the transport coefficient from the model of Braginskii (1965), and from Graille et al. (2009) respectively. Bold lines correspond to the case A, the other lines correspond to the case B.

In summary, we can conclude that the proposed method is verified for the fully ionized case. The main differences that are obtained between the two models are due to 1-the order of Laguerre Sonine polynomials that was used that is second order in Braginskii’s model, Balescu (1988), and third order in the proposed method, and 2- the nature of the collision operators used, Landau collision operators in the model of Braginskii, as opposed to Chapman and Cowling collision operators in the model of Graille et al. (2009). Additionally, the formulation of the transport properties that are considered in this paper are generalized for any type of partially ionized mixture.

5. Transport properties for a partially ionized helium-hydrogen plasma

5.1. Transport fluxes in thermo-chemical equilibrium

In order to simplify the analysis of the presented transport systems, we consider thermochemical equilibrium T_𝔢 = T_𝔥 = T, isobaric mixtures at rest. The total heat flux is entirely a function of the temperature gradient and magnetic field and may be written as

$$\begin{array}{c}\hfill {\mathit{q}}_{\mathfrak{h}}+{\mathit{q}}_{\mathfrak{e}}=-({\lambda }_{\mathfrak{h}}+{\overline{\overline{\lambda }}}_{\mathfrak{e}}+{\overline{\overline{\lambda }}}_S+{\overline{\overline{\lambda }}}_R){\mathit{\partial }}_{\mathit{x}}T,\end{array}$$(33)

where the Soret and reactive thermal conductivities may be written as

$$\begin{array}{c}\hfill {\overline{\overline{\lambda }}}_S=-p_{\mathfrak{e}}{\overline{\overline{\chi }}}_{\mathfrak{e}}{\theta }_{\mathfrak{e}}-{\displaystyle \sum _{j\in \mathrm{H}}[p_{\mathfrak{h}}{\chi }_{\mathfrak{h}j}+p_{\mathfrak{e}}{\overline{\overline{\chi }}}_{\mathfrak{e}}]{\theta }_i,}\end{array}$$(34)

$$\begin{array}{c}\hfill {\overline{\overline{\lambda }}}_R=-{\rho }_{\mathfrak{e}}{h}_{\mathfrak{e}}{\theta }_{\mathfrak{e}}-{\displaystyle \sum _{j\in \mathrm{H}}[{\rho }_j{h}_j+{\rho }_{\mathfrak{e}}{h}_{\mathfrak{e}}{\overline{\overline{\alpha }}}_{\mathfrak{e}j}]{\theta }_j,}\end{array}$$(35)

where θ_𝔢 and θ_i, i ∈ H are defined as

$$\begin{array}{c}\hfill {\theta }_{\mathfrak{e}}=-{\overline{\overline{D}}}_{\mathfrak{e}}[\frac{1}{x_{\mathfrak{e}}}\frac{\partial x_{\mathfrak{e}}}{\partial T}+\frac{{\overline{\overline{\chi }}}_{\mathfrak{e}}}{T}],\end{array}$$(36)

$$\begin{array}{c}\hfill {\theta }_i={\displaystyle \sum _{j\in \mathrm{H}}{D}_{\mathit{ij}}[\frac{1}{1-x_{\mathfrak{e}}}(\frac{\partial x_i}{\partial T}+\frac{\partial x_{\mathfrak{e}}}{\partial T}{\overline{\overline{\alpha }}}_{\mathfrak{e}j})+\frac{{\chi }_{\mathfrak{h}j}}{T}+\frac{p_{\mathfrak{e}}}{p_{\mathfrak{h}}}\frac{{\overline{\overline{\chi }}}_{\mathfrak{e}j}}{T}],i\in \mathrm{H}.}\end{array}$$(37)

Here, θ_𝔢 and θ_i correspond to diffusion velocities for a temperature gradient of 1, that are V_𝔢 = θ_𝔢∂_xT and V_i = θ_i∂_xT. We compute all the transport properties for the helium-hydrogen mixture S₁ for case A and case B.

Figures 8 and 9 present the parallel, perpendicular and transverse components of the electron thermal conductivity tensor ${\overline{\overline{\lambda }}}_{\mathfrak{e}}$ as a function of the temperature, for both cases. According to Fig. 8 (case A), the perpendicular component is equal to the parallel component for the entire range of temperatures, thus, the electron thermal conductivity is isotropic. Indeed, the pressure forces are dominating the magnetic pressure forces, so the plasma is unmagnetized. On the other hand, in Fig. 9 (case B), for temperatures higher than T = 5000 K, the electron thermal conductivity ${\overline{\overline{\lambda }}}_{\mathfrak{e}}$ is anisotropic since the magnitude of magnetic field is higher. This results in a transverse component that is higher than the perpendicular component of ${\overline{\overline{\lambda }}}_{\mathfrak{e}}$. Similar results have been obtained for the other anisotropic electron transport properties such as ${\overline{\overline{D}}}_{\mathfrak{e}}$ and ${\overline{\overline{\chi }}}_{\mathfrak{e}}$.

Fig. 8. Components of the electron thermal conductivity tensor ${\overline{\overline{\lambda }}}_{\mathfrak{e}}$ for the isotropic case A, at the third-order Laguerre Sonine polynomials, for the helium-hydrogen mixture S1 as a function of temperature.

Fig. 9. Components of the electron thermal conductivity tensor ${\overline{\overline{\lambda }}}_{\mathfrak{e}}$ for the anisotropic case B, at the third-order Laguerre Sonine polynomials, for the helium-hydrogen mixture S1 as a function of the temperature.

Figure 10 shows the heavy-particle thermal conductivity λ_𝔥, as a function of the temperature, for the case A and case B. In Fig. 10, strong differences between the two cases for a temperature higher than 6000 K can be seen. In the case A, λ_𝔥 increases from 1000 K to 9000 K, which is expected since λ_𝔥 is an increasing function of the temperature. However, in the case A after 9000 K, λ_𝔥 decreases. This decrease is due to the ionization of hydrogen. Indeed, the heavy particle thermal conductivity is related to a combination of the cross sections variations of all the heavy species in the mixture, which are proportional to the mole fractions of each heavy particles. This result is coherent with Fig. 1, which shows that the mole fraction of H is decreasing after 9000 K. Similar behavior as the case A have been observed for the case B, except that the ionization of H starts at 6000 K for this pressure. In Fig. 10, the second modulation observed around 12 000 K is due to the ionization of Helium as shown in Fig. 2.

Fig. 10. Heavy thermal conductivity λ𝔥, at the second-order of Laguerre Sonine approximation, for case A and case B, for the helium-hydrogen mixture S1 as a function of temperature.

Figure 11 shows the components of the total heat flux (Eq. (33)) as a function of the temperature, for the isotropic case A. It is clear that the reactive thermal conductivity ${\lambda }_R^{\parallel }$ is higher than the other components for certain ranges of temperature between 2200 K and 4300 K, where hydrogen dissociation occurs, and for temperature higher than 10 000 K, where hydrogen ionization takes place. The heavy thermal conductivity λ_𝔥 is the second term which dominates the total heat flux, and is higher than ${\lambda }_R^{\parallel }$ for a range of temperature between 4200 K and 10 000 K. The results here obtained are consistent with those of Scoggins et al. (2016).

Fig. 11. Component of the total heat flux (33) as a function of the temperature for the isotropic case A for the helium-hydrogen mixture S1.

Figure 12 and 13 show the parallel component of each term of the electron heavy-particle transport coefficients ${\alpha }_{\mathfrak{e}j}^{\Vert }$, as a function of the temperature, for the case A and case B. As before, each term of the electron heavy-particle transport tensor ${\alpha }_{\mathfrak{e}j}^{\Vert }$, is proportional to the mole fraction.

Fig. 12. Parallel component of the electron/heavy-particle transport coefficient ${\alpha }_{\mathfrak{e}j}^{\Vert },j\in \{\mathrm{He},{\mathrm{He}}^{+},\mathrm{H},{\mathrm{H}}_2,{\mathrm{He}}^{++},{\mathrm{H}}^{+}\}$, for the isotropic case A, for the helium-hydrogen mixture S1.

Fig. 13. Parallel component of the electron/heavy-particle transport coefficient ${\alpha }_{\mathfrak{e}j}^{\Vert },j\in \{\mathrm{He},{\mathrm{He}}^{+},\mathrm{H},{\mathrm{H}}_2,{\mathrm{He}}^{++},{\mathrm{H}}^{+}\}$, for the anisotropic case B, for the helium-hydrogen mixture S1.

5.2. Transport properties in a pore at the photosphere of the Sun

As done in the previous section, the transport coefficients of the previous helium-hydrogen mixture are computed for the conditions found in the upper layer of the solar convective zone from the radiative 3D MHD simulations of a pore by Kitiashvili et al. (2010). The simulation results are obtained for the computational domain of 6.4 × 6.4 × 5.5 Mm³ with the grid sizes: 50 × 50 × 43 km³, 25 × 25 × 21.7 km³ and 12.5 × 12.5 × 11 km³ (128² × 127, 256² × 253 and 512² × 505 mesh points). The domain includes an upper 5 Mm-deep layer of the convective zone and the solar chromosphere. In this section, the results and quantities are obtained from simulations of Kitiashvili et al. (2010), Wray et al. (2015, 2018) via a single-fluid model. This is a postprocessed calculation for quantities that belong to the multicomponent model presented. Since the electron heavy-particle collision frequency is high in the conditions chosen, a thermal equilibrium case T_𝔢 = T_𝔥 has been considered for the computation of the transport coefficients. For the sake of clarity, only results from a slice at a constant geometrical height z = −0.5 Mm, in the lower photosphere has been presented.

Figures 14–16 show snapshots of the distribution of the plasma beta parameter β_p, temperature T, and total mass density, respectively. As it can be seen, the temperature is varying from 4000 K to 6500 K, the plasma beta parameter is varying on a large range of magnitude, from weakly- to strongly-magnetized. In the snapshot of the simulation, a characteristic granulation pattern with the relatively hot (T >  5500 K) and less dense upflowing weakly-magnetized plasma in the middle of the granular cells can be observed. In addition, the lower temperature (T <  4500 K) and higher density downflowing strongly magnetized plasma at the intergranulation boundaries can be perceived (red lines of granulation). A strongly magnetized cold plasma can be seen in the middle of the snapshot.

Fig. 14. Plasma beta coefficient βp distribution from the radiative 3D MHD simulations of a pore by Kitiashvili et al. (2010).

Fig. 15. Temperature (K) distribution from the radiative 3D MHD simulations of a pore by Kitiashvili et al. (2010).

Fig. 16. Total mass density (kg m−3) distribution from the radiative 3D MHD simulations of a pore by Kitiashvili et al. (2010).

Figures 17–20 present the distribution of the heavy-particle heat flux λ_𝔥|∂_xT| and the ratios ${\lambda }_{\mathfrak{e}}^{\Vert }/{\lambda }_{\mathfrak{e}}^{\perp }$, ${\lambda }_R^{\parallel }$/λ_𝔥, and ${\lambda }_{\mathfrak{e}}^{\Vert }/{\lambda }_{\mathfrak{h}}$, respectively. Figure 18 shows that the electron thermal conductivity tensor ${\overline{\overline{\lambda }}}_{\mathfrak{e}}$ is almost isotropic everywhere, except in the middle of the snapshot where ${\lambda }_{\mathfrak{e}}^{\Vert }/{\lambda }_{\mathfrak{e}}^{\perp }=1.08$. In Fig. 19, the peak of reactive thermal conductivity is attributed to the dissociation reaction of hydrogen. Figure 20 shows that the electron thermal conductivity is small compared to the heavy thermal conductivity. These results are related to the results from Figs. 1 and 2 that show that the mole fraction of electrons is very small compared to the mole fraction of heavy particles in that range of temperature between 4000 K to 6500 K.

Fig. 17. Distribution of the heavy particle heat flux λ𝔥|∂xT|, for the helium-hydrogen mixture S1 based on the results of the radiative 3D MHD simulations of a pore by Kitiashvili et al. (2010).

Fig. 18. Ratio ${\lambda }_{\mathfrak{e}}^{\Vert }/{\lambda }_{\mathfrak{e}}^{\perp }$ distribution, computed at the third order of the Laguerre-Sonine polynomials approximation, for the helium-hydrogen mixture S1 based on the results of the radiative 3D MHD simulations of a pore by Kitiashvili et al. (2010).

Fig. 19. Ratio ${\lambda }_R^{\parallel }$/λ𝔥 distribution, computed at the third order of the Laguerre-Sonine polynomials approximation, for the helium-hydrogen mixture S1 based on the results of the radiative 3D MHD simulations of a pore by Kitiashvili et al. (2010).

Fig. 20. Ratio ${\lambda }_R^{\parallel }/{\lambda }_{\mathfrak{h}}$ distribution, computed at the third order of the Laguerre-Sonine polynomials approximation, for the helium-hydrogen mixture S2 based on the results of the radiative 3D MHD simulations of a pore by Kitiashvili et al. (2010).

5.3. Components of the generalized Ohm’s law in a pore at the Sun photosphere

As in Sect. 5.2, we compute the components of the generalized Ohm’s law from Eq. (25) using a helium-hydrogen mixture, from the simulation by Kitiashvili et al. (2010). According to the result found in Fig. 18, we assume an isotropic distribution of the transport properties. Figures 21–25 show the distribution of the resistive term, the electron battery term, the heavy-particle battery term and the Soret terms for electron and heavy particles respectively.

Fig. 21. Distribution of the resistive term (first term of Eq. (25)), computed at the third order of the Laguerre-Sonine polynomials approximation, for the helium-hydrogen mixture S1, based on the results of the radiative 3D MHD simulations of a pore by Kitiashvili et al. (2010).

Fig. 22. Distribution of the electron battery term (second term of Eq. (25)), computed at the third order of the Laguerre-Sonine polynomials approximation, for the helium-hydrogen mixture S1, based on the results of the radiative 3D MHD simulations of a pore by Kitiashvili et al. (2010).

Fig. 23. Distribution of the heavy battery term (third term of Eq. (25)), computed at the third order of the Laguerre-Sonine polynomials approximation, for the helium-hydrogen mixture S1, based on the results of the radiative 3D MHD simulations of a pore by Kitiashvili et al. (2010).

Fig. 24. Distribution of the Soret/Dufour term (fourth term of Eq. (25)), computed at the third order of the Laguerre-Sonine polynomials approximation, for the helium-hydrogen mixture S1, based on the results of the radiative 3D MHD simulations of a pore by Kitiashvili et al. (2010).

Fig. 25. Distribution of the Soret/Dufour term (last term of Eq. (25)), computed at the third order of the Laguerre-Sonine polynomials approximation, for the helium-hydrogen mixture S1, based on the results of the radiative 3D MHD simulations of a pore by Kitiashvili et al. (2010).

In the results obtained, the resistive term appears to be the higher term in the generalized Ohm’s law in both inside (magnitude 10²) and outside the pore (magnitude 1 − 10, not shown in Fig. 21). However, outside the pore, the battery term for heavy particles appears to be the second higher term. All the other terms are negligible inside and outside the pore (magnitude 10⁻⁷). In addition, the Soret and battery terms for electrons are negligible compared to the other terms of the generalized Ohm’s law. Indeed, this is due to the mole fraction of electrons which is very small compared to heavy particles under these conditions. These results are coherent with the mole fraction distribution presented in Figs. 1–3.

6. Conclusion

The present model, derived from the kinetic theory of Graille et al. (2009), is neither a single-fluid MHD nor a multi-fluid model. It is an intermediary model between the two that is a drift-diffusion multi-component model. Although the multi-fluid model and the multicomponent model are catching the same physics, their main difference is the scaling that has been used in the generalized Chapman-Enskog expansion. Their scaling leads to a thermal non-equilibrium multicomponent model with one momentum equation, where the electrons and each heavy particle diffuse in the hydrodynamic heavy-particle reference frame. These developments lead to an extended range of validity for partially and fully ionized plasma, as well as weakly, strongly and non-magnetized plasmas, and for a general multicomponent mixtures, which can be applied to conditions of the Sun’s atmophere. From the numerical point of view, the multicomponent model is prone to less stiffness compared to the conventional multi-fluid models, where source terms and coupling terms are considered in the momentum equations and energy equations of each particle. On the other hand, the single-fluid MHD model is shown to be a thermal equilibrium model with no relaxation term and decoupling between the particles. However, in this framework, this model shows less numerical stiffness in absence of source or relaxation terms. Fluid-models are valid only in cases where the plasma is in collisional regimes leading to small values of the Knudsen number. Thus, only perturbations of Maxwellian distribution functions are considered here.

From the set of governing equations, a generalized Ohm’s law has been derived. A general expression of the resistive term, as well as the battery term, has been obtained for a general multicomponent plasma. This general expression of the electric field can be simplified in a fully ionized plasma case (See Appendix B).

General conditions largely representative of the lower solar atmosphere have been chosen in order to compute all the transport properties for a helium-hydrogen mixture S₁. The latter are obtained by solving transport systems which are presented in Appendix C. A spectral Galerkin method based on a third-order Laguerre-Sonine polynomials approximation has been used. For convenience, we have implemented the transport model described in the paper in Mutation++, an open-source library¹. In order to validate our model, a comparison with the model of Braginskii (1965) has been performed in the case of a fully ionized plasma S₂. Both models are derived from the kinetic theory based on a Chapman-Enskog method. Differences are observed in 1-the structure of the governing equations and 2- the nature of the collision operators. While the heavy transport properties are anisotropic in Braginskii (1965), in our model, they remain isotropic. Nevertheless, under the chosen range of conditions, both behave as isotropic. In Braginskii (1965), the corresponding series of Laguerre-Sonine polynomials are truncated at the second order, whereas a third-order accuracy is reached in our model. Good agreement has been obtained for the considered fully ionized S₂ mixture in the chosen conditions.

Finally, using the Mutation++ library, the method allows us to compute all the transport properties for a partially ionized plasma for a given mixture. The obtained results strongly depend on the mole fraction between the species of the mixture. We have been able to identify the behavior of the transport coefficients related to the chemistry of the species in the partially ionized mixture S₁. Considering the complexity of the multicomponent model and the new terms that appear, postprocessing calculations were performed based on the results of a pore simulation for a partially ionized plasma of a mixture of Hydrogen and Helium S₁ in thermochemical equilibrium in the highly turbulent upper layer of the solar convective zone by Kitiashvili et al. (2010). These results allow us to both understand the effect and magnitude of these terms under the conditions we chose, especially the terms of the generalized Ohm’s law. Generally, the plasma being weakly ionized under these conditions, all the transport fluxes related to the electrons are negligible compared to those from heavy particles. In particular, we show a new formulation of the resistivity which depends on the multicomponent transport coefficients taking into account all possible interactions in the mixture. It appears to be much more general than standard resistivity such as ambipolar resistivity or Spitzer resistivity. Under conditions where the plasma is dominated by electromagnetic forces, which is the case in the center of the pore in Sect. 5, we find qualitatively that this term dominates the dynamics of the weakly ionized helium-hydrogen plasma. On the other hand, at the interboundary granulations (or more generally in the case of quiet sun photosphere conditions) when the plasma is weakly ionized and thermal pressure dominated, the new terms such as the Soret/Dufour effects (related to the temperature or pressure gradients) are dominating the plasma dynamics. In such conditions, these new terms allow magnetic reconnection process to occur. Indeed, according to Eq. (25), the parallel component of these terms may produce an electric field parallel to the magnetic field which may result in reconnection. These terms do not appear in the classical multi-fluid or single-fluid model (see Leake et al. 2013) since only a resistive term is generally considered, which allows magnetic reconnection to occur only in the magnetic pressure dominated plasma regime. Preliminary results have been obtained by Wargnier et al. (2018) for a fully ionized plasma in thermal nonequilibrium in the framework of a magnetic reconnection event.

Acknowledgments

The research of Q. Wargnier is funded by an Idex Paris-Saclay interdisciplinary IDI PhD grant, and relies on the the support of NASA Ames Research Center (ARC), Advanced Supercomputing Division, von Karman Institute for Fluid Dynamics, CMAP – Initiative HPC@Maths from Ecole Polytechnique and of Ecole doctorale de Mathématiques Hadamard. T. Magin is supported by a Jean d’Alembert chair of University Paris-Saclay. Part of this work was conducted during the 2018 NASA Summer Program at ARC. We would like to thank I.N. Kitiashvili for providing us with the data in order to evaluate the transport properties in Sect. 5.2. A. Alvarez Laguna and J. B. Scoggins were each funded by a post-doctoral fellowship from the Fondation Mathematique Jacques Hadamard (FMJH), LabEx Mathématique Hadamard - ANR11-LABX-0056.

References

Appendix A: Non-dimensional Boltzmann equations

In Graille et al. (2009), the non-dimensional Boltzmann equations for electron and heavy particles is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\partial }_t{\mathit{f}}_{\mathfrak{e}}& +\frac{1}{\epsilon M_{\mathfrak{h}}}({\mathit{C}}_{\mathfrak{e}}++\epsilon M_{\mathfrak{h}}{\mathit{v}}_{\mathfrak{h}})·{\mathit{\partial }}_{\mathit{x}}{\mathit{f}}_{\mathfrak{e}}\hfill \\ \hfill & +{\epsilon }^{-(1-b)}{\mathfrak{q}}_e({\mathit{C}}_{\mathfrak{e}}+\epsilon M_{\mathfrak{h}}{\mathit{v}}_{\mathfrak{h}})\otimes \mathit{B}·{\partial }_{{\mathit{C}}_{\mathfrak{e}}}{\mathit{f}}_{\mathfrak{e}}\hfill \\ \hfill & +(\frac{1}{\epsilon M_{\mathfrak{h}}}{\mathfrak{q}}_e\mathit{E}-\epsilon M_{\mathfrak{h}}\frac{D{\mathit{v}}_{\mathfrak{h}}}{\mathit{Dt}})·{\partial }_{{\mathit{C}}_{\mathfrak{e}}}{\mathit{f}}_{\mathfrak{e}}\hfill \\ \hfill & -({\partial }_{{\mathit{C}}_{\mathfrak{e}}}{\mathit{f}}_{\mathfrak{e}}\otimes {\mathit{C}}_{\mathfrak{e}}):{\mathit{\partial }}_{\mathit{x}}{\mathit{v}}_{\mathfrak{h}}=\frac{1}{{\epsilon }^2}{\mathcal{J}}_{\mathfrak{e}},\hfill \end{array}\end{array}$$(A.1)

and for each heavy species,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\partial }_t{\mathit{f}}_i& +\frac{1}{M_{\mathfrak{h}}}({\mathit{C}}_i+M_{\mathfrak{h}}{\mathit{v}}_{\mathfrak{h}})·{\mathit{\partial }}_{\mathit{x}}{\mathit{f}}_i\hfill \\ \hfill & +{\epsilon }^{-(1+b)}{\mathfrak{q}}_i({\mathit{C}}_i+\epsilon M_{\mathfrak{h}}{\mathit{v}}_{\mathfrak{h}})\otimes \mathit{B}·{\partial }_{{\mathit{C}}_i}{\mathit{f}}_i\hfill \\ \hfill & +(\frac{1}{\epsilon M_{\mathfrak{h}}}\frac{{\mathfrak{q}}_i}{m_i}\mathit{E}-M_{\mathfrak{h}}\frac{D{\mathit{v}}_{\mathfrak{h}}}{\mathit{Dt}}).{\partial }_{{\mathit{C}}_i}{\mathit{f}}_i\hfill \\ \hfill & -({\partial }_{{\mathit{C}}_i}{\mathit{f}}_i\otimes {\mathit{C}}_i):{\mathit{\partial }}_{\mathit{x}}{\mathit{v}}_{\mathfrak{h}}=\frac{1}{\epsilon }{\mathcal{J}}_i,i\in \mathrm{H}.\hfill \end{array}\end{array}$$(A.2)

Where the collision operators are defined as

$$\begin{array}{c}\hfill {\mathcal{J}}_{\mathfrak{e}}={\mathcal{J}}_{\mathfrak{e}\mathfrak{e}}({\mathit{f}}_{\mathfrak{e}},{\mathit{f}}_{\mathfrak{e}})+{\displaystyle \sum _{j\in \mathfrak{h}}{\mathcal{J}}_{\mathfrak{e}j}({\mathit{f}}_{\mathfrak{e}},{\mathit{f}}_j),}\end{array}$$(A.3)

$$\begin{array}{c}\hfill {\mathcal{J}}_i=\frac{1}{\epsilon }{\mathcal{J}}_{i\mathfrak{e}}({\mathit{f}}_i,{\mathit{f}}_{\mathfrak{e}})+{\displaystyle \sum _{j\in \mathfrak{h}}{\mathcal{J}}_{\mathit{ij}}({\mathit{f}}_i,{\mathit{f}}_j),i\in \mathrm{H}.}\end{array}$$(A.4)

Where C_i, i ∈ H and C_𝔢 are the peculiar velocities for heavy species and electron respectively.

Appendix B: Ohm’s law for a fully ionized plasma S₁

In a fully ionized plasma, the definition of the current density is

$$\begin{array}{c}\hfill \mathbf{I}=\mathfrak{n}\mathfrak{q}{\mathit{v}}_{\mathfrak{h}}+{\mathit{J}}_{\mathfrak{e}}\end{array}$$(B.1)

where the electron current density is defined by Eq. (4). It should be noted that in the particular case of a fully ionized plasma, no diffusion velocities of heavy species are considered, because the presented model is in the reference frame of heavy particles.

Then, similarly as the general case, the total current density can be expressed in function of fluxes, as follows,

$$\begin{array}{c}\hfill \mathbf{I}=\mathfrak{n}\mathfrak{q}{\mathit{v}}_{\mathfrak{h}}+\frac{{\left(n_{\mathfrak{e}}{\mathfrak{q}}_e\right)}^2}{p_{\mathfrak{e}}}{\overline{\overline{D}}}_{\mathfrak{e}}{\mathit{E}}^{\prime }-n_{\mathfrak{e}}{\mathfrak{q}}_e[{\overline{\overline{D}}}_{\mathfrak{e}}\frac{{\mathit{\partial }}_{\mathit{x}}{p}_{\mathfrak{e}}}{p_{\mathfrak{e}}}+{\overline{\overline{D}}}_{\mathfrak{e}}{\overline{\overline{\chi }}}_{\mathfrak{e}}{\mathit{\partial }}_{\mathit{x}}ln{T}_{\mathfrak{e}}]\end{array}$$(B.2)

It should be noted that the multicomponent electromagnetic matrices in a fully ionized plasma are

$$\begin{array}{cc}\hfill {\overline{\overline{M}}}_{{\mathit{E}}^{\prime }}& ={\overline{\overline{D}}}_{\mathfrak{e}},\hfill \end{array}$$(B.3)

$$\begin{array}{cc}\hfill {\overline{\overline{M}}}_{p_{\mathfrak{e}}}& ={\overline{\overline{D}}}_{\mathfrak{e}},\hfill \end{array}$$(B.4)

$$\begin{array}{cc}\hfill {\overline{\overline{M}}}_{p_j}& =0,j\in \mathrm{H}\hfill \end{array}$$(B.5)

$$\begin{array}{cc}\hfill {\overline{\overline{M}}}_{T_{\mathfrak{e}}}& ={\overline{\overline{D}}}_{\mathfrak{e}}{\overline{\overline{\chi }}}_{\mathfrak{e}}\hfill \end{array}$$(B.6)

$$\begin{array}{cc}\hfill {\overline{\overline{M}}}_{T_{\mathfrak{h}}}& =0.\hfill \end{array}$$(B.7)

Using the Maxwell-Ampere’s law ${\mathit{\partial }}_{\mathit{x}}\wedge \mathit{B}={\mu }_0\mathbf{I}$ in a non-relativistic context, the expression of the electric field as a function of fluxes can be obtained

$$\begin{array}{c}\hfill \mathit{E}=-{\mathit{v}}_{\mathfrak{h}}\wedge \mathit{B}+{\overline{\overline{\eta }}}_e{\mathit{J}}_{\mathfrak{e}}+\frac{{\mathit{\partial }}_{\mathit{x}}{p}_{\mathfrak{e}}}{{\mathfrak{n}}_{\mathfrak{e}}{\mathfrak{q}}_{\mathfrak{e}}}+\frac{p_{\mathfrak{e}}}{{\mathfrak{n}}_{\mathfrak{e}}{\mathfrak{q}}_{\mathfrak{e}}}{\overline{\overline{\chi }}}_{\mathfrak{e}}{\mathit{\partial }}_{\mathit{x}}ln{T}_{\mathfrak{e}}\end{array}$$(B.8)

where the electron resistivity tensor is ${\overline{\overline{\eta }}}_{\mathfrak{e}}={{\overline{\overline{D}}}_{\mathfrak{e}}}^{-1}{p}_{\mathfrak{e}}/{(n_{\mathfrak{e}}{q}_{\mathfrak{e}})}^2$.

Appendix C: Transport systems

C.1. Anisotropic transport systems

Anisotropic transport coefficients are computed in terms of solutions to complex linear systems. For μ ∈ {D_𝔢, (D_i)_i ∈ H}, transport systems read for Sonine polynomial order ξ ≥ 1

$$\begin{array}{cc}& {\displaystyle \sum _{q\in P_{\xi }}{L}_{\mathfrak{e}\mathfrak{e}}^{\mathit{pq}}{\alpha }_{\mathfrak{e}}^{q\mu (1)}(\xi )={\beta }_{\mathfrak{e}}^{p\mu },p\in P_{\xi },}\hfill \end{array}$$(C.1)

$$\begin{array}{cc}& {\displaystyle \sum _{q\in P_{\xi }}(L_{\mathfrak{e}\mathfrak{e}}^{\mathit{pq}}+i{L}_{\mathfrak{e}\mathfrak{e}}^{\mathit{Bpq}}){\alpha }_{\mathfrak{e}}^{q\mu (2)}(\xi )={\beta }_{\mathfrak{e}}^{p\mu },p\in P_{\xi },}\hfill \end{array}$$(C.2)

where P_ξ = {0, …, ξ − 1}. Likewise, for λ_𝔢, we have

$$\begin{array}{cc}& {\displaystyle \sum _{q\in P_{\xi 1}}{L}_{\mathfrak{e}\mathfrak{e}}^{\mathit{pq}}{\alpha }_{\mathfrak{e}}^{q{\lambda }_{\mathfrak{e}}(1)}(\xi )={\beta }_{\mathfrak{e}}^{p{\lambda }_{\mathfrak{e}}},p\in P_{\xi 1},}\hfill \end{array}$$(C.3)

$$\begin{array}{cc}& {\displaystyle \sum _{q\in P_{\xi 1}}(L_{\mathfrak{e}\mathfrak{e}}^{\mathit{pq}}+i{L}_{\mathfrak{e}\mathfrak{e}}^{\mathit{Bpq}}){\alpha }_{\mathfrak{e}}^{q{\lambda }_e(2)}(\xi )={\beta }_{\mathfrak{e}}^{p{\lambda }_{\mathfrak{e}}},p\in P_{\xi 1},}\hfill \end{array}$$(C.4)

where P_ξ1 = {1, …, ξ − 1} for ξ ≥ 2, the matrices $L_{\mathfrak{e}\mathfrak{e}}^{\mathit{pq}}$ and $L_{\mathfrak{e}\mathfrak{e}}^{\mathit{Bpq}}$ are given up to third order as

$$\begin{array}{cc}& L_{\mathfrak{e}\mathfrak{e}}^{\mathit{pq}}=\frac{16}{3}\frac{p}{k_B{T}_{\mathfrak{e}}}\sqrt{\frac{m_{\mathfrak{e}}}{2\pi k_B{T}_{\mathfrak{e}}}}{\stackrel{\sim }{L}}_{\mathfrak{e}\mathfrak{e}}^{\mathit{pq}},\hfill \end{array}$$(C.5)

$$\begin{array}{cc}& {\stackrel{\sim }{L}}_{\mathfrak{e}\mathfrak{e}}^{00}={\displaystyle \sum _{j\in \mathrm{H}}{x}_j{\overline{Q}}_{\mathfrak{e}j}^{(1,1)},}\hfill \end{array}$$(C.6)

$$\begin{array}{cc}& {\stackrel{\sim }{L}}_{\mathfrak{e}\mathfrak{e}}^{01}={\stackrel{\sim }{L}}_{\mathfrak{e}\mathfrak{e}}^{10}={\displaystyle \sum _{j\in \mathrm{H}}{x}_j(\frac{5}{2}{\overline{Q}}_{\mathfrak{e}j}^{(1,1)}-3{\overline{Q}}_{\mathfrak{e}j}^{(1,2)}),}\hfill \end{array}$$(C.7)

$$\begin{array}{cc}& {\stackrel{\sim }{L}}_{\mathfrak{e}\mathfrak{e}}^{02}={\stackrel{\sim }{L}}_{\mathfrak{e}\mathfrak{e}}^{20}={\displaystyle \sum _{j\in \mathrm{H}}{x}_j(\frac{35}{8}{\overline{Q}}_{\mathfrak{e}j}^{(1,1)}-\frac{21}{2}{\overline{Q}}_{\mathfrak{e}j}^{(1,2)}+6{\overline{Q}}_{\mathfrak{e}j}^{(1,3)}),}\hfill \end{array}$$(C.8)

$$\begin{array}{cc}& {\stackrel{\sim }{L}}_{\mathfrak{e}\mathfrak{e}}^{11}={\displaystyle \sum _{j\in \mathrm{H}}{x}_j(\frac{25}{4}{\overline{Q}}_{\mathfrak{e}j}^{(1,1)}-15{\overline{Q}}_{\mathfrak{e}j}^{(1,2)}+12{\overline{Q}}_{\mathfrak{e}j}^{(1,3)})+x_{\mathfrak{e}}\sqrt{2}{\overline{Q}}_{\mathfrak{e}\mathfrak{e}}^{(2,2)},}\hfill \end{array}$$(C.9)

$$\begin{array}{cc}\hfill {\stackrel{\sim }{L}}_{\mathfrak{e}\mathfrak{e}}^{12}& ={\stackrel{\sim }{L}}_{\mathfrak{e}\mathfrak{e}}^{21}={\displaystyle \sum _{j\in \mathrm{H}}{x}_j(\frac{175}{16}{\overline{Q}}_{\mathfrak{e}j}^{(1,1)}-\frac{315}{8}{\overline{Q}}_{\mathfrak{e}j}^{(1,2)}+57{\overline{Q}}_{\mathfrak{e}j}^{(1,3)}-30{\overline{Q}}_{\mathfrak{e}j}^{(1,4)}})\hfill \\ \hfill & +x_{\mathfrak{e}}\sqrt{2}(\frac{7}{4}{\overline{Q}}_{\mathfrak{e}\mathfrak{e}}^{(2,2)}-2{\overline{Q}}_{\mathfrak{e}\mathfrak{e}}^{(2,3)}),\hfill \end{array}$$(C.10)