© The Authors 2025

1. Introduction

Owing to the relatively cool temperatures in the solar photosphere and chromosphere, the plasma in these layers (primarily consisting of hydrogen) is partially ionised. Moving into the upper layers, the magnetic pressure exceeds the gas pressure and the magnetic field determines the behaviour of charges. Neutral particles can only sense the presence of the magnetic field through collisions with charges. The friction resulting from charge-neutral collisions causes the dissipation of kinetic and magnetic energies, which are then converted into internal energy, raising the plasma temperature. This process becomes more effective when both the magnetohydrodynamic (MHD) and collision timescales are comparable (Zaqarashvili et al. 2011a; Soler et al. 2013a,b; Ballester et al. 2018; Popescu Braileanu & Keppens 2021). Furthermore, Khodachenko et al. (2004, 2006) demonstrated that, under conditions observed in prominences, wave damping resulting from charge-neutral collisions is more pronounced than damping caused by viscosity or thermal conductivity. In addition, Shelyag et al. (2016) found that charge-neutral collisional damping is more efficient than Ohmic dissipation for propagating waves in flux tubes.

Depending on the scales involved, describing the interaction between charges and neutrals might require a model that includes both components as separated fluids. One significant reference with respect to the description of a multi-fluid plasma (N-fluid, or NF) is the seminal paper of Braginskii (Braginskii 1965), among other works (see e.g. Chapman & Cowling 1970; Schunk & Nagy 2009). This approach requires us to computie the moments of the Boltzmann equation for the distribution function of each component of the plasma. The Chapman-Enskog successive-approximation technique calculates the high-order moments of the distribution (Schunk 1977; Meier & Shumlak 2012). Braginskii’s method is applicable to highly collisional plasmas and is valid across most of the chromosphere. For plasmas with intermediate or low collision coupling, a more general model such as the one developed by Hunana et al. (2019) and Hunana et al. (2022) is required. Recently, Hunana (2025) extended Braginskii’s set of equations to a 22-moment model for a 2F plasma composed by neutrals and charges.

When dealing with partially ionised solar plasmas, a common method for simulating the interaction between charged and neutral particles is using a single-fluid model (1F; e.g. Zweibel & Brandenburg 1997; Ballester et al. 2018; Nóbrega-Siverio et al. 2020a). This model includes ambipolar diffusion among other non-ideal terms in the generalised Ohm’s law. Outside of solar physics, it has been used, for example, in simulating neutron stars (Jones 1987), the interstellar medium (Brandenburg 2019), or weakly ionised protoplanetary disks (Bai & Stone 2011). Many numerical studies in solar physics have invoked the ambipolar term. For instance, Cheung & Cameron (2012), Martínez-Sykora et al. (2012) and Khomenko et al. (2018) found ambipolar diffusion to produce measurable effects in magneto-convection simulations, exceeding artificial diffusivities at chromospheric layers. Martínez-Sykora et al. (2016) analysed the misalignment of magnetic field lines with the visible structure of spicules and due to the interaction between charges and neutrals, invalidating the MHD assumption that plasma motion consistently follows the same magnetic field lines. Nóbrega-Siverio et al. (2020b) and Martínez-Sykora et al. (2020a, 2023) found that nonequilibrium (NEQ) conditions may affect the efficiency of ambipolar diffusion. Recently, Gómez-Míguez et al. (2024) incorporated inelastic collisions into ambipolar diffusion term definition, finding the corrections to the ambipolar coefficient to be negligible in quiet Sun conditions. Works by Khodachenko et al. (2004, 2006), Forteza et al. (2006), Song & Vasyliūnas (2011), Soler et al. (2015a), Shelyag et al. (2016), Cally & Khomenko (2019), Khomenko & Cally (2019), and Popescu Braileanu et al. (2021), among others, have investigated the damping of MHD waves and currents due to ambipolar diffusion. They found high-frequency Alfvén and fast waves to be effectively damped through this mechanism in many different scenarios such as quiet Sun, flux tubes, or prominences. However, for the case of slow waves, Forteza et al. (2006) found them to be weakly affected by ion-neutral collisions. Moreover, Terradas et al. (2002, 2005), Carbonell et al. (2006), Forteza et al. (2008), and Soler et al. (2015a) found them to be more affected by other non-ideal mechanisms, such as viscosity or thermal conduction.

In recent years, many NF studies have appeared in the context of Solar physics. In terms of wave theory, Zaqarashvili et al. (2011a, 2013) and Soler et al. (2013a,b) characterised the 2F normal modes in an uniform atmosphere, focussing on the cutoff regions not present in the 1F approach (see e.g. Kulsrud & Pearce 1969; Kamaya & Nishi 1998). Expanding the analysis of 2F waves, Cally (2023) studied the 10 modes resulting from the linear equations applied to a homogeneous atmosphere, focussing on the efficiency in generating those modes and Cally & Gómez-Míguez (2023) applied a ray tracing formalism for the study of fast-to-slow mode conversion in 2F under the assumption of weak dissipation. We refer to the review by Ballester et al. (2018) and Soler (2024) for recent developments in 2F waves.

Regarding numerical simulations, Popescu Braileanu et al. (2019a,b) or Zhang et al. (2021) studied wave heating in a 1D stratified atmosphere for different wave periods, finding the temperature increments to decrease for shorter wave periods. Snow & Hillier (2021) delved into 1D shock propagation, finding neutrals to stabilise the wavefront for finite collisions. Furthermore, Meier (2011), Maneva et al. (2017), Ballai (2019), and Zhang et al. (2021) expanded this research by incorporating ionisation-recombination effects, revealing challenges in configuring a static background atmosphere when hydrostatic equilibrium is imposed. Recently, Snow et al. (2023) also included radiative ionisation-recombination rates for a (5 + 1 levels) hydrogen atom for investigating shock propagation and Popescu Braileanu & Keppens (2024) added radiative cooling under coronal conditions to study the growth of the thermal instability and acoustic wave damping. Apart from 2F modeling, Martínez-Gómez et al. (2017) explored a 5F approach that includes hydrogen and helium, while Martínez-Gómez et al. (2018) used a 3F model to study the non-linear regime of Alfvén waves in a 1D homogeneous atmosphere. In addition, the Ebysus code (Martínez-Sykora et al. 2020b) has been shown to be capable of performing NF simulations for different species, treating ionisation and excitation states separately.

Although NF models offer greater flexibility, make fewer assumptions, and thus offer a more accurate description of chromospheric plasma at small scales, it is still possible that a simpler 1F model is sufficiently precise for application in the chromosphere, depending on the frequency of the studied process. Hence, it is necessary to compare both approaches and verify the applicability of the 1F approach under conditions similar to those of the solar chromosphere. This study aims to enhance our understanding of the heating mechanisms associated with interactions between charged particles and neutrals and to determine the specific scenarios where the 1F approximation becomes less accurate than the more detailed 2F approximation. For this purpose, we employed the MANCHA-2F code developed by Popescu Braileanu et al. (2019a) to solve both the 1F and 2F equations and compare their solutions in a simple scenario involving the vertical propagation of fast magneto-acoustic waves in a solar-like chromospheric model, following the methodology of previous research (e.g. Popescu Braileanu et al. 2019b).

The paper is organised as follows. Section 2 outlines the mathematical models utilised in this study. Section 3 presents the results of the numerical experiments on the propagation of fast magneto-acoustic waves. Section 4 presents the analytical solution for the wave equation. Section 5 analyses the reasons for the 1F and 2F discrepancies. Finally, Section 6 summarises our findings.

2. Physical description of the plasma

In this work, we assume a hydrogen plasma consisting of electrons, protons, and neutral particles, without any other species, molecules, or negative hydrogen ions. We focus only on the effects due to elastic charge-neutral collisions. We do not include inelastic collisions such as ionisation or recombination and we do not consider viscosity or heat flux. By restricting ourselves to a 5-moment model (Schunk 1977), we are able to consistently apply these assumptions in a closed model. We note that Ohmic heating, the Hall effect, or Biermann battery effect are not included in the generalised Ohm’s law.

2.1. The two-fluid (2F) model

Our model follows a 2F approach, where elastic collisions produce a coupling between a neutral fluid with a charge fluid. The neutral fluid has only one component, while the charge fluid consists of strongly coupled protons and electrons. In this instance, we are examining a 2F model comprising charges and neutrals, as established by Meier & Shumlak (2012), Leake et al. (2013), or Popescu Braileanu (2020). The set of moment equations is

$$\begin{array}{cc}\hfill \frac{\partial {\rho }_{\mathrm{n}}}{\partial t}+\mathrm{\nabla }·({\rho }_{\mathrm{n}}{\mathbf{u}}_{\mathrm{n}})& =0,\hfill \end{array}$$(1a)

$$\begin{array}{cc}\hfill \frac{\partial {\rho }_{\mathrm{c}}}{\partial t}+\mathrm{\nabla }·({\rho }_{\mathrm{c}}{\mathbf{u}}_{\mathrm{c}})& =0,\hfill \end{array}$$(1b)

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}({\rho }_{\mathrm{n}}{\mathbf{u}}_{\mathrm{n}})+\mathrm{\nabla }·({\rho }_{\mathrm{n}}{\mathbf{u}}_{\mathrm{n}}\otimes {\mathbf{u}}_{\mathrm{n}}+{\widehat{p}}_{\mathrm{n}})& ={\rho }_{\mathrm{n}}\mathbf{g}+{\mathbf{R}}_{\mathrm{n}},\hfill \end{array}$$(1c)

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}({\rho }_{\mathrm{c}}{\mathbf{u}}_{\mathrm{c}})+\mathrm{\nabla }·({\rho }_{\mathrm{c}}{\mathbf{u}}_{\mathrm{c}}\otimes {\mathbf{u}}_{\mathrm{c}}+{\widehat{p}}_{\mathrm{c}}+{\widehat{p}}_{\mathrm{m}})& ={\rho }_{\mathrm{c}}\mathbf{g}-{\mathbf{R}}_{\mathrm{n}},\hfill \end{array}$$(1d)

$$\begin{array}{cc}\hfill \frac{\partial e_{\mathrm{n}}}{\partial t}+\mathrm{\nabla }·[e_{\mathrm{n}}{\mathbf{u}}_{\mathrm{n}}+{\widehat{p}}_{\mathrm{n}}·{\mathbf{u}}_{\mathrm{n}}]& ={\rho }_{\mathrm{n}}\mathbf{g}·{\mathbf{u}}_{\mathrm{n}}+M_{\mathrm{n}},\hfill \end{array}$$(1e)

$$\begin{array}{cc}\hfill \frac{\partial e_{\mathrm{c}}}{\partial t}+\mathrm{\nabla }·[e_{\mathrm{c}}{\mathbf{u}}_{\mathrm{c}}+({\widehat{p}}_{\mathrm{c}}+{\widehat{p}}_{\mathrm{m}})·{\mathbf{u}}_{\mathrm{c}}]& ={\rho }_{\mathrm{c}}\mathbf{g}·{\mathbf{u}}_{\mathrm{c}}-M_{\mathrm{n}},\hfill \end{array}$$(1f)

where the usual thermodynamic variables ${\rho }_{\alpha },{\mathbf{u}}_{\alpha },{\widehat{p}}_{\alpha },and{e}_{\alpha }$ are, respectively, the mass density, the centre of mass velocity, the stress tensor, and the total energy of the fluid α = c, n for charges and neutrals, respectively. We note that the charges include the sum of protons and electrons. The scalar product is denoted by ‘⋅’ and the outer product by ‘⊗’. The total energy is, in each case,

$$\begin{array}{c}\hfill e_{\mathrm{n}}=e_{\mathrm{int},\mathrm{n}}+\frac{1}{2}{\rho }_{\mathrm{n}}{u}_{\mathrm{n}}^2;e_{\mathrm{c}}=e_{\mathrm{int},\mathrm{c}}+\frac{1}{2}{\rho }_{\mathrm{c}}{u}_{\mathrm{c}}^2+p_{\mathrm{m}};e_{\mathrm{int},\alpha }=\frac{p_{\alpha }}{\gamma -1},\end{array}$$(2)

with e_int, α the internal energy and the adiabatic index γ = 5/3. Density and pressure obey the ideal gas relations

$$\begin{array}{cc}\hfill p_{\mathrm{n}}& =\frac{{\rho }_{\mathrm{n}}}{m_H}{k}_{\mathrm{B}}{T}_{\mathrm{n}},\hfill \end{array}$$(3a)

$$\begin{array}{cc}\hfill p_{\mathrm{c}}& =\frac{2{\rho }_{\mathrm{c}}}{m_{\mathrm{H}}}{k}_{\mathrm{B}}{T}_{\mathrm{c}},\hfill \end{array}$$(3b)

where $p_{\alpha }=\frac{1}{3}\mathrm{Tr}{\widehat{p}}_{\alpha }$ is the scalar pressure, m_H the mass of hydrogen atom, k_B the Boltzmann constant, and T_α the temperature.

The gravitational acceleration, g, and the magnetic field, B, characterise the macroscopic forces acting on the system. The magnetic force is described throughout the divergence of the magnetic stress tensor, ${\widehat{p}}_m$, defined as

$$\begin{array}{cc}\hfill {\widehat{p}}_m& =p_m\widehat{I}-\frac{1}{{\mu }_0}\mathbf{B}\otimes \mathbf{B},\hfill \end{array}$$(4a)

$$\begin{array}{cc}\hfill p_m& =\frac{B^2}{2{\mu }_0},\hfill \end{array}$$(4b)

where μ₀ is the vacuum magnetic permeability, p_m is the magnetic pressure, and $\widehat{I}$ is the 3 × 3 identity tensor. The magnetic force is commonly expressed in terms of the current density, $\mathbf{J}=\frac{1}{{\mu }_0}\mathbf{\nabla }\times \mathbf{B}$, so $\mathrm{\nabla }·{\widehat{p}}_{\mathrm{m}}=-\mathbf{J}\times \mathbf{B}$.

The magnetic field evolves according to the induction equation

$$\begin{array}{c}\hfill \frac{\partial \mathbf{B}}{\partial t}+\mathrm{\nabla }·({\mathbf{u}}_{\mathrm{c}}\otimes \mathbf{B}-\mathbf{B}\otimes {\mathbf{u}}_{\mathrm{c}})=0,\end{array}$$(5)

and obeys the no monopole condition

$$\begin{array}{c}\hfill \mathrm{\nabla }·\mathbf{B}=0.\end{array}$$(6)

The terms R_n and M_n describe the elastic collisional interaction between the fluids (Braginskii 1965). The elastic collision term for the mass equation is null since no particles are produced or destroyed in elastic processes. For the case of a hydrogen plasma, the collisional terms are given by

$$\begin{array}{cc}\hfill {\mathbf{R}}_{\mathrm{n}}& =\alpha {\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}({\mathbf{u}}_{\mathrm{c}}-{\mathbf{u}}_{\mathrm{n}}),\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill M_{\mathrm{n}}& =\frac{1}{2}\alpha {\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}(u_{\mathrm{c}}^2-u_{\mathrm{n}}^2)+\frac{1}{\gamma -1}\frac{k_{\mathrm{B}}}{m_{\mathrm{H}}}\alpha {\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}(T_{\mathrm{c}}-T_{\mathrm{n}}),\hfill \end{array}$$(8)

where α is the collision parameter,

$$\begin{array}{c}\hfill \alpha (T_{\mathrm{eff}})=\frac{1}{m_{\mathrm{H}}^2}\sqrt{\frac{8k_{\mathrm{B}}{T}_{\mathrm{eff}}}{\pi }}(\sqrt{m_{\mathrm{pn}}}{\sigma }_{\mathrm{pn}}+\sqrt{m_{\mathrm{en}}}{\sigma }_{\mathrm{en}}).\end{array}$$(9)

The subindex ‘p’ refers to protons, while m_βn = m_βm_n/(m_β + m_n), T_eff = (T_c + T_n)/2, and σ_βn = 10⁻¹⁹ m² are the cross-sections, assuming hard sphere elastic collisions (Draine et al. 1983; Leake et al. 2012), for β = p, e. Using these definitions, the collision frequencies are

$$\begin{array}{cc}\hfill {\nu }_{\mathit{cn}}& =\alpha {\rho }_n,\hfill \end{array}$$(10a)

$$\begin{array}{cc}\hfill {\nu }_{\mathit{nc}}& =\alpha {\rho }_c.\hfill \end{array}$$(10b)

These frequencies define the collision time scales τ_cn = 1/ν_cn and τ_nc = 1/ν_nc. If the MHD time scale of the studied processes is similar to or smaller than the collision ones, the 2F nature of the plasma arises.

2.2. A single fluid (1F) model with ambipolar diffusion

Following the same notation as in the 2F Section 2.1, the 1F set of equations is expressed as

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·(\rho \mathbf{u})& =0,\hfill \end{array}$$(11a)

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}(\rho \mathbf{u})+\mathrm{\nabla }·(\rho \mathbf{u}\otimes \mathbf{u}+\widehat{p}+{\widehat{p}}_{\mathrm{m}})& =\rho \mathbf{g},\hfill \end{array}$$(11b)

$$\begin{array}{cc}\hfill \frac{\partial e}{\partial t}+\mathrm{\nabla }·[e\mathbf{u}+(\widehat{p}+{\widehat{p}}_{\mathrm{m}})·\mathbf{u}+{\mathbf{S}}^{\ast }]& =\rho \mathbf{g}·\mathbf{u},\hfill \end{array}$$(11c)

$$\begin{array}{cc}\hfill \frac{\partial \mathbf{B}}{\partial t}+\mathrm{\nabla }·(\mathbf{u}\otimes \mathbf{B}-\mathbf{B}\otimes \mathbf{u})& =-\mathrm{\nabla }\times {\mathbf{E}}^{\ast },\hfill \end{array}$$(11d)

with the thermodynamical relations,

$$\begin{array}{cc}\hfill e& =e_{\mathrm{int}}+\frac{1}{2}\rho u^2+p_{\mathrm{m}};e_{\mathrm{int}}=\frac{p}{\gamma -1}\hfill \end{array}$$(12a)

$$\begin{array}{cc}\hfill p& =\frac{\rho }{\stackrel{\sim }{\mu }{m}_{\mathrm{n}}}{k}_{\mathrm{B}}T,\hfill \end{array}$$(12b)

$$\begin{array}{cc}\hfill \stackrel{\sim }{\mu }& =\frac{1}{2-{\xi }_{\mathrm{n}}},\hfill \end{array}$$(12c)

$$\begin{array}{cc}\hfill {\xi }_{\mathrm{c}}+{\xi }_{\mathrm{n}}& =\frac{{\rho }_{\mathrm{c}}}{\rho }+\frac{{\rho }_{\mathrm{n}}}{\rho }=1.\hfill \end{array}$$(12d)

The amount of neutrals and charges is encapsulated on the mass fraction, ξ_α, and on the average particle mass, $\stackrel{\sim }{\mu }$. In the 1F approach, ambipolar diffusion is included through the generalised Ohm’s law. This effect is present in Eqs. (11c) and (11d) by introducing the non-ideal electric field E^* and the Poynting flux S^*, which can be expressed as

$$\begin{array}{cc}\hfill {\mathbf{E}}^{\ast }& =-{\eta }_{\mathrm{A}}\frac{(\mathbf{J}\times \mathbf{B})\times \mathbf{B}}{B^2}=-{\stackrel{\sim }{\eta }}_{\mathrm{A}}\mathbf{B}\times \mathrm{\nabla }·{\widehat{p}}_{\mathrm{m}},\hfill \end{array}$$(13a)

$$\begin{array}{cc}\hfill {\mathbf{S}}^{\ast }& =\frac{1}{{\mu }_0}{\mathbf{E}}^{\ast }\times \mathbf{B}=-\frac{{\eta }_{\mathrm{A}}}{{\mu }_0}\frac{(\mathbf{J}\times \mathbf{B})\times \mathbf{B}}{B^2}\times \mathbf{B}\hfill \end{array}$$(13b)

$$\begin{array}{cc}& =-\frac{1}{{\mu }_0}({\stackrel{\sim }{\eta }}_{\mathrm{A}}\mathbf{B}\times \mathrm{\nabla }·{\widehat{p}}_{\mathrm{m}})\times \mathbf{B},\hfill \end{array}$$(13c)

where we express the ambipolar coefficient as

$$\begin{array}{c}\hfill {\stackrel{\sim }{\eta }}_{\mathrm{A}}=\frac{{\xi }_{\mathrm{n}}^2}{\alpha {\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}}·\end{array}$$(14)

We also introduced the more usual ambipolar coefficient ${\eta }_{\mathrm{A}}=B^2{\stackrel{\sim }{\eta }}_{\mathrm{A}}$ and expressed the ambipolar terms in a more standard fashion. However, we made use of the non-standard notation, as explained in the following sections.

The momentum exchange terms in the 2F (see Eq. (7)) and the ambipolar diffusion refer to the same process: the diffusion of neutrals and charges. It is worthwhile to clarify the connection between these two concepts. Here, we follow the steps outlined by, namely, Khomenko et al. (2014).

We consider a 1F model that accounts for velocity drifts between neutrals and charges,

$$\begin{array}{c}\hfill \mathbf{w}={\mathbf{u}}_{\mathrm{c}}-{\mathbf{u}}_{\mathrm{n}}.\end{array}$$(15)

In this approach, the dynamics of the plasma is described by the centre of mass of the components, which moves at the velocity

$$\begin{array}{c}\hfill \mathbf{u}=\frac{{\rho }_{\mathrm{n}}{\mathbf{u}}_{\mathrm{n}}+{\rho }_{\mathrm{c}}{\mathbf{u}}_{\mathrm{c}}}{{\rho }_{\mathrm{n}}+{\rho }_{\mathrm{c}}}·\end{array}$$(16)

With the above considerations, we can relate the drift velocity of neutrals and charges respect to the centre of mass, w_α, with the drift velocity w, by using Eqs. (15) and (16),

$$\begin{array}{cc}\hfill {\mathbf{w}}_{\mathrm{n}}& ={\mathbf{u}}_{\mathrm{n}}-\mathbf{u}=-{\xi }_{\mathrm{c}}\mathbf{w},\hfill \end{array}$$(17a)

$$\begin{array}{cc}\hfill {\mathbf{w}}_{\mathrm{c}}& ={\mathbf{u}}_{\mathrm{c}}-\mathbf{u}={\xi }_{\mathrm{n}}\mathbf{w}.\hfill \end{array}$$(17b)

From the inspection of the charge energy equation and the 2F induction equation (see Eqs. (1f) and (5)), it follows that the magnetic field terms are related to the velocity of the charged fluid. The non-ideal terms E^* and S^* (see Eqs. (13a) and (13b)) result from the change of the charge velocity to the velocity of the centre of mass frame velocity, so

$$\begin{array}{c}\hfill {\mathbf{u}}_{\mathrm{c}}=\mathbf{u}+{\xi }_{\mathrm{n}}\mathbf{w},\end{array}$$(18)

and,

$$\begin{array}{c}\hfill {\mathbf{E}}^{\ast }=\mathbf{E}+\mathbf{u}\times \mathbf{B}=-{\xi }_{\mathrm{n}}\mathbf{w}\times \mathbf{B}.\end{array}$$(19)

The drift velocity, w, is not a native variable of the 1F model, as it requires us to evolve the velocities for neutrals and charges separately. Therefore, it is necessary to make some assumptions in order to account for the drift velocity. By combining the momentum equations for neutrals and charges (Eqs. (1c) and (1d)), we find

$$\begin{array}{cc}& {\xi }_{\mathrm{n}}(\frac{\partial }{\partial t}({\rho }_{\mathrm{c}}{\mathbf{u}}_{\mathrm{c}})+\mathrm{\nabla }·({\rho }_{\mathrm{c}}{\mathbf{u}}_{\mathrm{c}}\otimes {\mathbf{u}}_{\mathrm{c}}))\hfill \\ \hfill & -{\xi }_{\mathrm{c}}(\frac{\partial }{\partial t}({\rho }_{\mathrm{n}}{\mathbf{u}}_{\mathrm{n}})+\mathrm{\nabla }·({\rho }_{\mathrm{n}}{\mathbf{u}}_{\mathrm{n}}\otimes {\mathbf{u}}_{\mathrm{n}}))=-{\xi }_{\mathrm{n}}\mathrm{\nabla }·{\widehat{p}}_{\mathrm{m}}-\mathbf{G}-\alpha {\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}\mathbf{w},\hfill \end{array}$$(20)

similar to Eq. (129) by Khomenko et al. (2014). The terms on the left hand side (LHS) are the inertial terms, which is approximately expressed as

$$\begin{array}{cc}& {\xi }_{\mathrm{n}}(\frac{\partial }{\partial t}({\rho }_{\mathrm{c}}{\mathbf{u}}_{\mathrm{c}})+\mathrm{\nabla }·({\rho }_{\mathrm{c}}{\mathbf{u}}_{\mathrm{c}}\otimes {\mathbf{u}}_{\mathrm{c}}))\hfill \\ \hfill & -{\xi }_{\mathrm{c}}(\frac{\partial }{\partial t}({\rho }_{\mathrm{n}}{\mathbf{u}}_{\mathrm{n}})+\mathrm{\nabla }·({\rho }_{\mathrm{n}}{\mathbf{u}}_{\mathrm{n}}\otimes {\mathbf{u}}_{\mathrm{n}}))\approx \rho {\xi }_{\mathrm{c}}{\xi }_{\mathrm{n}}\frac{\partial \mathbf{w}}{\partial t}·\hfill \end{array}$$(21)

The ratio between the inertial terms and the collisional momentum exchange term R_n can be roughly estimated by

$$\begin{array}{c}\hfill {\varpi }_{\mathrm{inertia}}=\frac{\omega }{\alpha \rho },\end{array}$$(22)

with ω as the characteristic frequency of the process.

If ϖ_inertia ≪ 1, the inertial terms can be neglected in Eq. (20), then

$$\begin{array}{c}\hfill {\mathbf{w}}_{1\mathrm{F}}=\frac{{\stackrel{\sim }{\eta }}_{\mathrm{A}}}{{\xi }_{\mathrm{n}}^2}({\xi }_{\mathrm{n}}\mathbf{J}\times \mathbf{B}-\mathbf{G})=-\frac{{\stackrel{\sim }{\eta }}_{\mathrm{A}}}{{\xi }_{\mathrm{n}}^2}({\xi }_{\mathrm{n}}\mathrm{\nabla }·{\widehat{p}}_{\mathrm{m}}+\mathbf{G}).\end{array}$$(23)

The first term on the right-hand side (RHS) of Eq. (23) is the drift velocity due to magnetic forces, which is usually referred to as ambipolar diffusion (see also the discussion by Nóbrega-Siverio et al. 2020a). This contribution is included in the 1F model presented here through the terms E^* and S^* (see Eqs. (13a) and (13b)). The second one is the thermal pressure function, given by

$$\begin{array}{c}\hfill \mathbf{G}={\xi }_{\mathrm{n}}\mathrm{\nabla }·{\widehat{p}}_{\mathrm{c}}-{\xi }_{\mathrm{c}}\mathrm{\nabla }·{\widehat{p}}_{\mathrm{n}},\end{array}$$(24)

which describes the drift due to the thermal pressure forces. In the absence of viscosity and assuming thermal equilibrium (T_n = T_c) and scalar pressures, it can be reduced to

$$\begin{array}{c}\hfill \mathbf{G}={\xi }_{\mathrm{n}}(1-{\xi }_{\mathrm{n}})\mathrm{\nabla }(\stackrel{\sim }{\mu }p)-(1+{\xi }_{\mathrm{n}})\stackrel{\sim }{\mu }p\mathrm{\nabla }{\xi }_{\mathrm{n}},\end{array}$$(25)

which is equivalent to the expression obtained by Forteza et al. (2006). In principle, we ignore this term in the 1F approach, since we study a weakly ionised atmosphere and G → 0 in both the limits of weak and strong ionisation (see the discussion by Khodachenko et al. 2004, 2006; Pinto et al. 2008). Additionally, in Eq. (23), there is a missing term proportional to the current density and to the electron mass (see e.g. Braginskii 1965; Khomenko et al. 2014; Benavides & Flierl 2020). This term comes from the drift velocity between ions and electrons and it is neglected in the 2F model discussed in this work (for more details, see the model decription by Popescu Braileanu 2020).

2.3. Internal energy equation

The energy equation can be cast in several alternative forms, showing explicitly different physical properties of the system. The internal energy equation is specially useful because it is connected with the temperature throughout the state equation. In Equations (1e), (1f) and (11c), it is possible to remove the explicit contribution of the kinetic and magnetic energy by using the momentum and induction equations, obtaining the following for the 2F model:

$$\begin{array}{cc}\hfill \frac{\partial e_{\mathrm{int},\mathrm{n}}}{\partial t}+\mathrm{\nabla }·(e_{\mathrm{int},\mathrm{n}}{\mathbf{u}}_{\mathrm{n}})& =-p_{\mathrm{n}}\mathrm{\nabla }·{\mathbf{u}}_{\mathrm{n}}+Q_{\mathrm{n}},\hfill \end{array}$$(26a)

$$\begin{array}{cc}\hfill \frac{\partial e_{\mathrm{int},\mathrm{c}}}{\partial t}+\mathrm{\nabla }·(e_{\mathrm{int},\mathrm{c}}{\mathbf{u}}_{\mathrm{c}})& =-p_{\mathrm{c}}\mathrm{\nabla }·{\mathbf{u}}_{\mathrm{c}}+Q_{\mathrm{c}}.\hfill \end{array}$$(26b)

Hereafter, we simplified the notation by using only the scalar pressure. The RHS terms represent the internal energy sources. The first one is the adiabatic source term and Q_α accounts for the energy variation due to elastic collisions. Explicitly, they are

$$\begin{array}{cc}\hfill Q_{\mathrm{n}}& =\frac{1}{2}\alpha {\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}{w}^2+\frac{1}{\gamma -1}\frac{k_{\mathrm{B}}}{m_{\mathrm{H}}}\alpha {\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}(T_{\mathrm{c}}-T_{\mathrm{n}}),\hfill \end{array}$$(27a)

$$\begin{array}{cc}\hfill Q_{\mathrm{c}}& =\frac{1}{2}\alpha {\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}{w}^2-\frac{1}{\gamma -1}\frac{k_{\mathrm{B}}}{m_{\mathrm{H}}}\alpha {\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}(T_{\mathrm{c}}-T_{\mathrm{n}}).\hfill \end{array}$$(27b)

The first term in Eqs. (27a) and (27b) is called the frictional heating (FH) (see e.g. Draine et al. 1983; Draine 1986; Leake et al. 2014; Ni et al. 2018). It depends on the squared modulus of the drift velocity between charges and neutrals. The second term is the thermal exchange (TE). While the FH is always positive leading to the same heating rate for both components, TE sign can vary and it is the opposite for neutrals and charges (implying that it tends to homogenise the temperature of both species). The sum of Q_n and Q_c

$$\begin{array}{c}\hfill Q_{\mathrm{n}}+Q_{\mathrm{c}}=\alpha {\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}{w}^2\end{array}$$(28)

accounts for the variation of the total internal energy due to the collisions, finding that it is strictly positive, as the TE terms cancel each other. Then, we can say that differences in velocity between neutrals and charges contribute to increase the total internal energy of the plasma and, therefore, its temperature, but temperature differences between charges and neutrals do not (see e.g. the discussion by Popescu Braileanu 2020).

Regarding the 1F equations, from Eqs. (11b) and (11c), it follows that

$$\begin{array}{cc}\hfill \frac{\partial e_{\mathrm{int}}}{\partial t}+\mathrm{\nabla }·(e_{\mathrm{int}}\mathbf{u})& =-p\mathrm{\nabla }·\mathbf{u}+Q,\hfill \end{array}$$(29a)

$$\begin{array}{cc}\hfill Q& ={\eta }_{\mathrm{A}}{J}_{\perp }^2={\stackrel{\sim }{\eta }}_{\mathrm{A}}{(\mathrm{\nabla }·{\widehat{p}}_{\mathrm{m}})}^2,\hfill \end{array}$$(29b)

where ${\mathbf{J}}_{\perp }=\mathbf{J}·(\widehat{I}-\mathbf{b}\otimes \mathbf{b})$ is the component of the current density perpendicular to the magnetic field, with $\mathbf{b}=\frac{\mathbf{B}}{B}$. Again, we find that Q is strictly positive, indicating that the internal energy of the centre of mass will increase, similar to the total internal energy in the 2F approach.

Fig. 1. Representation of the characteristic parameters of the atmosphere. In the top panel, the density of neutrals is in solid black line and the density of charges is in dashed-dotted line. Their scale is shown on the left axis. The right axis labels represent the neutral-charge (solid red line) and the charge-neutral (dashed-dotted red line) collision frequencies. The characteristic speeds are displayed in the bottom panel: the neutral (solid black line), charge (dashed line) and 1F (dashed-dotted green line) sound speeds and the 1F (solid yellow line) and charge (‘dashed-dot’ yellow line) Alfvén speeds. The vertical dashed line indicates the equipartition layer βplasma = 1.

3. Vertical propagation of MA waves in a stratified atmosphere with a homogeneous horizontal magnetic field. Numerical experiments

In this work, we solve both 2F and 1F equations using the code MANCHA-2F (Popescu Braileanu et al. 2019a) based on the code MANCHA3D (Khomenko & Collados 2006; Modestov et al. 2024). Collision times vary by orders of magnitude through the chromosphere and in the deep layers they are much smaller than the typical time scales of hydrodynamic processes. Short collision times require small integration time steps for an explicit numerical scheme to satisfy the CFL condition. To handle this issue, MANCHA-2F introduces a semi-implicit scheme for updating the contribution of the collision terms in the 2F model, speeding up the computational time of the simulations. The code can also optionally solve the linearised system of equations, offering the possibility of separately assessing the linear and non-linear effects.

The numerical experiments are performed in a stratified atmosphere along the vertical axis (z-axis). The base is located at the bottom of the chromosphere, about 500 km above the solar surface. In this region, hydrogen is the primary donor of free electrons, which allows for a relatively self-consistent evaluation of charge density in our model. The top boundary is at the base of the transition region, where temperature starts to grow strongly.

It must be noted that the approximations underlying the 2F model require that the electron-proton collisional frequency is above the proton-neutral collisional frequency, the condition which is only fulfilled, on average, starting from the upper photosphere or the low chromosphere upwards (see Figure 1 in Khomenko et al. 2014), as well as Vranjes & Krstic (2013), Soler et al. (2015b), Ballester et al. (2018). In our work we consider heights starting from 600 km upwards, where this condition is mostly verified. Additionally, for the low wave frequencies considered in our work, it will be shown that both 1F and 2F models practically converge to the ideal MHD at the bottom part of the simulation domain (Piddington 1956). The application of both 1F and 2F models throughout the whole atmosphere is needed for consistency of their comparison. The ion-electron decoupling effects are not included in the 1F model either, so this aspect is not a source of discrepancies for the present discussion.

We set the equipartition region (where β_plasma = p/p_m = 1) at the middle of the physical domain. The background magnetic field is homogeneous and is oriented along the x-axis. We take the value of the background magnetic field to be ∼17 G to set the equipartition region at the middle of the domain. The hydrostatic equilibrium is computed from a VAL3C temperature profile (Vernazza et al. 1981), taking the mass density of hydrogen and electrons at the base of the chromosphere. This configuration resembles the one from Popescu Braileanu et al. (2019b), except that the magnetic field is considered constant with height in the present work. Using the model with horizontal magnetic field allows us to reduce the problem to 1D propagation while still considering mode conversion, since the waves will change their nature from acoustic to magnetic on their way upwards. In the Sun, horizontal magnetic fields are frequently found in the canopy regions surrounding intergranular or supergranular magnetic field concentrations (Wiegelmann et al. 2014).

We imposed a hydrostatic equilibrium (HS), in which each fluid is balanced by its own forces. As the background magnetic field is homogeneous, it does not contribute to the force balance. Then, pressure profiles are given by

$$\begin{array}{cc}\hfill p_{\alpha ,0}(z)& =p_{\alpha ,0}(0)exp(-{\displaystyle {\int }_0^z}\frac{dz\prime }{H_{\alpha }(z\prime )}),\hfill \end{array}$$(30a)

$$\begin{array}{cc}\hfill H_{\mathrm{n}}& =\frac{k_{\mathrm{B}}{T}_{\mathrm{n},0}}{m_{\mathrm{H}}g};H_{\mathrm{c}}=\frac{2k_{\mathrm{B}}{T}_{\mathrm{c},0}}{m_{\mathrm{H}}g};H=\frac{k_{\mathrm{B}}{T}_0}{{\stackrel{\sim }{\mu }}_0{m}_{\mathrm{H}}g},\hfill \end{array}$$(30b)

where H_α is the pressure scale height and the label ‘0’ refers to the model atmosphere. We omit the index α for the 1F variables. For a uniform magnetic field, w₀ = 0, making the equilibrium states equivalent for both 2F and 1F models. However, if we consider the case of a non force-free magnetic field, i.e. $\mathrm{\nabla }·{\widehat{p}}_{\mathrm{m}}\ne 0$, then w₀ ≠ 0 (see discussion by Gómez-Míguez et al. 2024).

The characteristic speeds of the atmosphere are the sound speeds c_sα and the Alfvén speeds c_Ac and c_A, given by

$$\begin{array}{c}\hfill c_{\mathrm{s}}=\sqrt{\frac{\gamma p_0}{{\rho }_0}};c_{\mathrm{s}\alpha }=\sqrt{\frac{\gamma p_{\alpha ,0}}{{\rho }_{\alpha ,0}}};c_{\mathrm{A}}=\frac{B_{\mathrm{x},0}}{\sqrt{{\mu }_0{\rho }_0}};c_{\mathrm{Ac}}=\frac{B_{\mathrm{x},0}}{\sqrt{{\mu }_0{\rho }_{\mathrm{c},0}}}.\end{array}$$(31)

The speeds c_s and c_A are representative for the highly coupled plasma, whereas c_sα and c_Ac apply to the neutral and charged fluids when they are non-coupled. Here, we focus on the first case.

Fig. 1 illustrates the atmospheric variables. The collision frequencies mainly depend on the density, decreasing exponentially with height. The hydrostatic equilibrium determines the ionisation fraction, which takes into account only the initial conditions at the base and does not include ionisation/recombination effects. The ionisation fraction (which coincides with the charge mass fraction) varies with height due to the different pressure scales of neutrals and charges from 10⁻⁴ at the bottom of our atmosphere to 10⁻² at its top. Due to the low ionisation fraction in this model atmosphere, 1F profiles of thermodynamic variables and characteristic speeds closely resemble those for neutrals. The sound speeds are nearly constant, reflecting the fact that the temperature varies smoothly. The Alfvén speeds exhibit the most significant differences due to the low number of charges. However, if the components of the plasma are highly coupled, only the 1F Alfvén and sound speeds are relevant for the study of waves (Piddington 1956; de Pontieu & Haerendel 1998; Kumar & Roberts 2003; Soler et al. 2013a,b).

For the numerical experiments, we launched a vertically propagating fast magneto-acoustic wave of a certain period from the bottom of the atmosphere. With MANCHA-2F, we solve both linear and non-linear systems of equations. Unless specified otherwise, we are referring to the non-linear solutions hereafter. The driver introduces the pure fast magneto-acoustic mode as a perturbation by fixing the amplitude of u_c, z to be 10⁻³c_s, taking the local value of the sound speed at the bottom layer. The other variables are then evaluated according to the analytical solution. Details regarding the driver configuration can be found in Appendix A. At the bottom, charges and neutrals are highly coupled, which means they move at the same velocity, u_z (Zaqarashvili et al. 2011a; Soler et al. 2013b). Therefore, we set the amplitude of u_z in the 1F driver to be the same as that in the 2F driver. We performed a set of 2F and 1F simulations, varying the wave period and we select the cases of periods P = 1, 3, and 5 s, whose frequency is close to the collision frequency at the high layers of the atmosphere (about 2 s⁻¹).

The understanding of the numerical effects in the simulations is key for this analysis, especially since 1F and 2F employ different numerical schemes (explicit vs semi-implicit). In favour of performing a fair comparison, we restrict ourselves to the study of low-amplitude waves, despite the amount of heating we produce is negligible when compared to the background temperature. Low amplitudes allow for low non-linear effects, and we do not need to include any hyper-diffusivity or shock diffusivity for controlling numerical effects. We fix the time step to be the same for 1F and 2F simulations, choosing in all the cases the most restrictive one. We also perform tests for different spatial resolutions, finding that for the cases of 5 s and 3 s periods, dz ≈ 500 m is enough, but the 1 s period requires 4 times finer resolution. For the numerical stability reasons, a 6th order low-pass filter is applied to the solution each certain amount of simulation time, which is the same for 1F and 2F simulations. In all the experiments, we set a Perfectly Matched Layer (PML) (Berenger 1994; Parchevsky & Kosovichev 2007) of the same size and parameters, at the top boundary. In addition to these numerical configurations, we also study the effects of the numerical scheme in the solutions. The 1F code uses a 4th order Runge-Kutta method for the integration, while the 2F code uses a semi-implicit method for the collisional terms, which is at most second order (for more details, see Popescu Braileanu et al. 2019a). The 2F code has the option for computing the solution using an explicit scheme, the same as the 1F code. We have performed 2F simulations with both schemes, finding no significant differences between the solutions for our resolutions. Based on all of these arguments, we are confident that the differences between the solutions are caused by physical differences between the models.

The results of the simulations are shown in Fig. 2. At the bottom of the atmosphere, the gas pressure dominates over magnetic forces, and the waves behave essentially like acoustic-gravity waves. Since the waves move from more dense layers to less dense layers, the amplitude of velocity perturbation grows along the wave path. At a certain height, wave dissipation due to charge-neutral collisions becomes sufficiently effective in overcoming growth due to stratification and the amplitude of the waves starts to decay. As noted by Popescu Braileanu et al. (2019a,b) or Zhang et al. (2021), we also observe that damping reduces the non-linear effects of the wavefront. As expected, collisional damping is more effective for short-period waves, since their frequencies are closer to the collision frequencies.

Fig. 2. Groups of panels from top to bottom: comparison between the numerical solution of velocity, temperature and magnetic field perturbations for waves of period P = 1 s, 3 s and 5 s at 250 s after the start of the simulation. Solid red lines: Perturbations in the velocity, temperature of charges and perturbations in the magnetic field in the 2F model. Solid cyan lines: Same for neutrals in 2F model. Black dashed lines: Same for the 1F model. The equipartition layer βplasma = 1 is indicated with a grey dashed vertical line.

Broadly speaking, the 1F fluid solution agrees with the 2F results. As shown by Popescu Braileanu et al. (2019b), the 1F solution reproduces essentially the neutral fluid’s behaviour, since neutrals dominate in number. In line with their findings, we also observe a change in phase between charges and neutrals at the upper layers, leading to a velocity drift between the species. This difference is significant for the velocities, while temperatures remain similar for charges and neutrals because of the TE. Magnetic field oscillations are also very close for both models.

When analysing the solutions depending on the wave period, it leads to the conclusion that long-period waves exhibit greater decoupling effects, which may seem surprising. More in detail, we represent in Fig. 3 the centre of mass velocity and the drift velocity from the 2F simulations, computed as Eqs. (16) and (15). We find u ≫ w in most of the domain, except for the layers where the wave is more damped, finding there that u ≈ w. Due to the strong wave damping, it is less evident for a 1 s wave period. However, a zoom-in inspection of the oscillations above the equipartition region reveals u ≈ w. As stated before, the wave of larger period is the one leading to the largest drift velocity.

Fig. 3. Comparison of the centre of mass velocity (u, cyan curves) and drift velocity (w, magenta curves) for the simulations represented in Fig. 2, at 250 s of simulated time. The values were computed from the 2F simulations. The vertical dashed line indicates the equipartition layer βplasma = 1.

The 1F model is a valid description of plasma phenomena with characteristic timescales exceeding those of collisions. This model assumes that plasma components rapidly adjust their momentum and energy, behaving as a single fluid. When the timescales become closer to, but are still larger than collision times, ambipolar diffusion and other non-ideal effects due to neutrals may become relevant to describe the plasma behaviour in the 1F model. This limit was explored for 2F waves by Zaqarashvili et al. (2011a) and Soler et al. (2013a,b), who performed a normal mode analysis for a homogeneous atmosphere, and more recently Popescu Braileanu et al. (2019b) extended the study to a vertically stratified atmosphere. Our numerical experiments show that increasing the wave period in the range of periods we have explored leads to larger values of the drift velocity at the upper layers. This was already shown by other authors as Popescu Braileanu et al. (2019b), where they compared the standard deviation of w_z for different wave periods. This result can be clarified by noting that shorter-period waves are more effectively damped, preventing them from reaching the upper layers where the most notable variations between the two models are expected. Conversely, large period waves can reach upper layers, where the fluids are less coupled. Therefore, it is not only important to check for the ratio ω/ν_nc when studying decoupling effects, but also consider the efficiency of the wave for transporting energy to layers of low coupling. If we continue increasing the wave period, the trend must change because both the 2F and 1F models must converge to the ideal MHD model (see e.g. Zaqarashvili et al. 2011a; Soler et al. 2013a,b; Cally & Gómez-Míguez 2023). We do not explore this question in this work, but it can be anticipated from Fig. 6 by Popescu Braileanu et al. (2019b), who studied the propagation of waves of periods from 1 s up to 20 s, finding that waves about 5 s period are the ones producing the largest drift velocities.

Several works, namely, Braginskii (1965), Zaqarashvili et al. (2011a), or Khomenko et al. (2014) have derived the set of 1F equations from multifluids, obtaining terms proportional to w_α. Equations (17a) and (17b) show that for the case of a pure hydrogen plasma, w_α ∝ w (it also stands if there are more components and all the neutrals move at velocity u_n and all the charges at velocity u_c). We also have that for our weakly ionised model atmosphere, w_n → 0 and w_c ≈ w (see Eqs. (17a) and (17b)). Fig. 3 shows that u and w become comparable as the waves reach upper layers, and discrepancies between 1F and 2F models could appear. We will discuss this point in more detail in Section 5.

One of our primary objectives is to compare the heating effects in the 2F and 1F models. In this context, we represent in Fig. 4 (left column) the source terms resulting from collisional and ambipolar heating of each model, as defined by Eqs. (27a), (27b), and (29b), respectively. We multiply them by the factor

$$\begin{array}{c}\hfill N_{\mathrm{T}}=\frac{(\gamma -1)m_{\mathrm{H}}\stackrel{\sim }{\mu }}{\rho k_{\mathrm{B}}},\end{array}$$(32)

following the work by Popescu Braileanu et al. (2019b), so they are expressed in Ks⁻¹. At the bottom of the domain, the plasma is highly coupled and the heating terms play little role, and the wave behaves as expected from the ideal MHD model. As the waves propagate upwards, the collisional damping becomes more and more significant and the amplitude of the heating terms rises up until they reach a maximum. The peak position is higher for longer periods, consistently with the view of Osterbrock (1961). When comparing these terms with the wave behaviour shown in Fig. 2, we observe that the position of the heating peaks are located after the velocity amplitude starts to be damped. The wave amplitude continues decaying after the peak, and the wave has less energy to release. Consequently, the increment of temperature decays as well. Most of the differences between the 2F and 1F solutions are observed in the upper layers for longer-period waves (see left panels in Fig. 4). For the 1 s period wave, we observe a transition in the behavior of heating rates. In deep layers, 2F heating has a higher amplitude than the 1F one and moving up, the 1F heating is dominating. Ambipolar heating is driven by magnetic forces, and the good agreement between model results proves that ambipolar diffusion is the main contributor to the drift velocity, hence, the impact of the G term (see Eq. (23)) on the solution for the atmosphere tested here is negligible. Furthermore, magnetic forces are able to drive larger drift velocities between neutrals and charges than pressure forces, as consequence of the low ionisation fraction of the atmosphere (Khodachenko et al. 2004, 2006; Pinto et al. 2008).

Fig. 4. Left column: Sum of the 2F collisional heating terms (solid red line, see Eqs. (27b) and (27a)) and the 1F ambipolar heating term (black dashed line, see Eq. (29b))), at 190 s of simulated time. The factor NT (see Eq. (32)) is introduced in all the cases. Right column: Time-average temperature for the 1F and 2F models (solid and dashed lines, respectively). Different time intervals are represented with different colours. The vertical axis range is the same in all the panels. The 2F curves are computed from Eq. (34). From top to bottom, the rows refer to the 1 s, 3 s, and 5 s wave periods. The averaging time intervals are the same for the cases of 1 s and 5 s. We also include residuals (1F–2F) to better compare the differences between the models. The vertical dashed line indicates the equipartition layer βplasma = 1.

Next, we focussed on the key variable, the temperature. We can compute the average temperature increase by averaging the temperature perturbation over a wave period for each layer of the atmosphere. For the case of the 1-second wave period, we chose to average several periods to reduce the noise. We search for the peaks and valleys of the wave and construct the wave envelopes by interpolation. Then, we compute the mean value of the top and bottom envelopes and, finally, the time average. For the 2F wave, we compute the centre of mass temperature, so temperatures of both models can be compared directly. By combining Eqs. (12b) and (B.1c), it follows

$$\begin{array}{c}\hfill T=\frac{\stackrel{\sim }{\mu }{m}_{\mathrm{H}}}{\rho k_{\mathrm{B}}}{\displaystyle \sum _{\{n,c\}}}(p_{\alpha }+\frac{1}{3}{\rho }_{\alpha }{w}_{\alpha }^2)=\frac{\stackrel{\sim }{\mu }{m}_{\mathrm{H}}}{\rho k_{\mathrm{B}}}{\displaystyle \sum _{\{n,c\}}}{p}_{\alpha }+\frac{1}{3}\frac{\stackrel{\sim }{\mu }{m}_{\mathrm{H}}}{k_{\mathrm{B}}}{\xi }_{\mathrm{c}}{\xi }_{\mathrm{n}}{w}^2,\end{array}$$(33)

where we apply Eqs. (17a) and (17b) on the last step.

We can safely neglect the second term in Eq. (33) because ξ_c ≪ 1 and ρw² ≪ ∑p_α, being null at the equilibrium. Then, we find

$$\begin{array}{c}\hfill T\approx \stackrel{\sim }{\mu }({\xi }_{\mathrm{n}}{T}_{\mathrm{n}}+2{\xi }_{\mathrm{c}}{T}_{\mathrm{c}}).\end{array}$$(34)

This expression is similar to the one used, for example, by Popescu Braileanu et al. (2019b). Those converge for a low ionisation fraction, but may lead to different values when the plasma is highly ionised.

The comparison of the centre of mass temperatures from each model is shown in the right panels of Fig. 4, for different time ranges. It is apparent that the average temperature increases over time in both models. For all the wave periods, differences between the 2F and 1F models are minor, but they become more evident at heights above the peak in temperature increase. Below the peak, differences are subtle but still discernible. After the peak, the 1F model gives a more significant average temperature increase than the 2F model. We observe the opposite behaviour at heights before reaching the T₁ peak. The same pattern applies for all the considered periods. We have excluded the upper layers from the calculation for the 1-second period wave because the wave amplitude is too small and the computation becomes dominated by noise. It is not straightforward to compute the temperature rate per time for each layer from the heating terms and also the method for computing the average temperature may introduce some systematics. However, we can state that both the heating terms and the average temperature increase exhibit a similar behaviour but, quantitatively, we cannot confirm if the computations match.

Having examined how the temperature changes as a function of height in the atmosphere, we study next how the average temperature increase changes as a function of time. In Fig. 5, we present the average temperature increase at the heating peak (different for the three periods) over time for each of the simulations. The short-period wave data points overlap, indicating a good agreement between the 2F and 1F models. Conversely, there are more differences visible for the longer-period waves. Since the wave energy is proportional to the wave amplitude squared, if there is more energy available, then it is possible to increase more the internal energy of the plasma, even if the dissipation is less effective. These results are consistent with Fig. 14 by Zhang et al. (2021), where they represented ⟨T₁⟩ for different wave periods, but they were not including a magnetic field in their study.

Fig. 5. Evolution of the average temperature increase at the heating peak as a function of time. Diamonds: P = 5 s; dots: P = 3 s; stars: P = 1 s. The 2F data points are represented in black and the 1F ones in red.

It is important to highlight the linear dependence of ⟨T₁⟩ on time, in agreement with what was previously reported by Popescu Braileanu et al. (2019b) from their simulations with a similar set-up. This linear dependence is consequence of launching waves of small amplitude, meaning that energy dissipation is solely dependent on the initial background state and the wave frequency (Braginskii 1965; Priest 2014; Tracy et al. 2014). This increment is minor compared to the background temperature, however, one has to take into account that we use very low initial amplitudes, about 10 m s⁻¹. In the solar case, at the upper photosphere, the wave amplitudes reach ∼1000 m s⁻¹ (see, for example, Beck et al. 2009), so ⟨T₁⟩ will be proportionally larger.

The analysis of the heating terms represented in Fig. 4 brings another curious feature. We observe two groups of peaks whose envelopes have different amplitudes and shapes. This behaviour is more clear when looking at the case of 1 s wave period. From Eqs. (27a) and (27b), it follows that Q_n + Q_c ∝ w_z². We show in Fig. 6 the shape of w_z² for all the wave periods simulated, obtaining that w_z² exhibits the same behaviour as Q. We have also computed w_z² for the same set-up, but solving the linear equations with MANCHA-2F. The comparison shows that the two groups are present in both cases, but in the linear simulation they follow the same envelope. Back to Fig. 3, if we look at one of the peaks of u_z, the slope of the wave is steeper after the peak than before, for those more sharped. As consequence, the peaks and valleys in w_z become asymmetric due to the non-linear effects. Once the amplitude of w_z² starts to decay, the non-linear effects become weaker and the amplitude of the two groups of peaks tend to the linear amplitude. It is interesting to observe that for the case of 1 s wave period, the envelopes of w_z² have the maximum values at different layers. A similar result was obtained by Popescu Braileanu et al. (2019b) (see Fig. 8) for the same set-up. It is also possible to guess the same behaviour for the 3 s case. Despite the small amplitude of the waves, w_z is highly sensitive to the non-linear effects due to their high frequency.

Fig. 6. Comparison of wz2 for a linear (cyan line) and non-linear (magenta line) 2F simulations, at 250 s of simulated time. From top to bottom, the wave periods are P = 1, 3 and 5 s, respectively. The vertical dashed line indicates the equipartition layer βplasma = 1.

We completed our analysis of the simulated waves by looking into the energy distribution. From Braginskii (1965), Bray & Loughhead (1974). Walker (2004, 2014), Cally & Gómez-Míguez (2023) or Cally (2023), the wave energy terms can be defined as

$$\begin{array}{c}\hfill W_{\mathrm{K},\alpha }=\frac{1}{2}{\rho }_{\alpha ,0}{u}_{\alpha }^2;W_{\mathrm{ac},\alpha }=\frac{p_{\alpha ,1}^2}{2{\rho }_{\alpha ,0}{c}_{\mathrm{s}\alpha }^2};W_{\mathrm{m}}=\frac{B_1^2}{2{\mu }_0}·\end{array}$$(35)

From left to right, these expressions refer to the kinetic energy density, the acoustic energy density, and the magnetic energy density. Compared with the works cited above, we do not include the buoyancy term here as we have ω ≫ ν_ac, where the acoustic cutoff frequency reaches a maximum value of ${\nu }_{\mathrm{ac}}=\frac{c_{\mathrm{s}}}{2H}\approx 5$ mHz at the bottom of our atmosphere, which is located around the temperature minimum. Therefore, we consider total energy W_T to be only the sum of the kinetic, acoustic and magnetic energy, for all the fluid components considered in each model. Figs. 7 and 8 show the amount of energy of each type that the waves are carrying along the model atmosphere. Those terms are multiplied by the 1F group velocity of the fast mode, $v_{\mathrm{g}}=\sqrt{c_{\mathrm{s}}^2+c_{\mathrm{A}}^2}$, so the total energy times the group velocity (the total wave energy flux) is constant in the 1F limit in absence of dissipation mechanisms, based on the wave-action conservation (see e.g. Tracy et al. 2014; Walker 2014). All the waves exhibit a constant energy flux at lower layers of the atmosphere, where the neutrals and charges are highly coupled by collisions. Once the collisional dissipation becomes effective, the total energy decreases. This was also noted by Zhang et al. (2021) when analysing the acoustic energy of the waves. We also find that the distance from the layer where the wave starts to decay up to the layer where the wave effectively vanishes decreases for short wave periods, which evidences the more efficiency of charge-neutral collisions in dissipation high-frequency waves.

Fig. 7. Representation of the 1F wave energy flux contributions for the three wave periods studied in this work. From the top to the bottom, the wave periods are 1 s, 3 s, and 5 s, respectively. Blue line: Magnetic energy (Wm). Green line: Acoustic energy (Wac). Red line: Kinetic energy (WK). Black line: Total energy (WT). The vertical dashed line indicates the equipartition layer βplasma = 1.

Fig. 8. Representation of the 2F wave energy flux contributions for the three wave periods studied in this work. From the top to the bottom in groups of two plots, the wave periods are 1 s, 3 s, and 5 s, respectively. Blue line: the magnetic energy (Wm). Green line: Acoustic energy (Wac, α). Red line: Kinetic energy (WK, α). Black line: Total energy (WT). The terms of the neutrals are shown in solid lines and terms for the charges in dashed lines. We split the terms WK, c and Wac, c in panels separated from neutrals due to the difference in scale. The vertical dashed line indicates the equipartition layer βplasma = 1.

When looking at the behaviour of the energy components in 1F (see Fig. 7), we observe that, at the lower layers, the wave energy is divided into kinetic and acoustic energy almost in the same proportion, with a small contribution of the magnetic energy. This behavior is the print of an acoustic wave, whose energy equally splits into acoustic and kinetic (see e.g. Landau & Lifshitz 1987; Walker 2014). As the wave propagates upwards, the kinetic energy keeps constant but the acoustic energy gradually exchanges to magnetic energy. At the point where energy dissipation becomes efficient, all the energies start to decay, but the exchange between acoustic to magnetic energy keeps the same trend. In the neighborhood of the equipartition region, the wave reaches a layer where acoustic and magnetic energies are in the same proportion, and as the wave continues its path, the acoustic energy is gradually exchanged into magnetic energy. In other words, the analysis of the wave energy shows that the wave transforms from acoustic to magnetic (see e.g. Schunker & Cally 2006; Khomenko & Collados 2006; Walker 2014; Cally & Gómez-Míguez 2023).

If we follow the same analysis for the 2F waves (see Fig. 8), we find that the energy terms related to the neutrals and magnetic energy show approximately the same behaviour as in the 1F waves. Regarding the kinetic and acoustic energy of charges, we find them to be orders of magnitude smaller than those of the neutrals because of the low ionisation fraction. In addition, at deep layers, those energy fluxes are clearly different; the acoustic energy approximately doubles the kinetic energy. This can be explained as follows. In those layers, the ideal MHD solution is still valid and the phase speed is essentially c_s, therefore, ω ≈ k_zc_s. Then, from Eq. (A.3), it follows that p_c, 1 ≈ ρ_c, 0c_su_c, z, where we have neglected the growth term since we are only interested in a small region and also assume that u_c, z = u_z, in light of Fig. 3. Introducing this relation into the expression for the acoustic energy for charges in (35), we obtain

$$\begin{array}{c}\hfill W_{\mathrm{ac},\mathrm{c}}\approx \frac{p_{\mathrm{c},0}}{2p_0}{\rho }_0{u}_{c,z}^2\approx {\rho }_{c,0}{u}_{c,z}^2,\end{array}$$(36)

which is the double of the kinetic energy of the charges. In the last step, we have applied the limit of weakly ionisation. Unlike the case of the neutrals, W_K, c and W_ac, c grow up at the deep layers and then decay after peaking, much similar to the behavior of W_m.

For all the wave periods considered in this work, we observe that the dissipation becomes more efficient as more energy is exchanged into magnetic energy. Forteza et al. (2006) found the fast waves to be efficiently damped by charge-neutral collisions in prominences. In their study, the Alfvén speeds were above 10 km s⁻¹, while the temperature was 8000 K. For this temperature in the VAL3C model, the sound speed is less than 10 km s⁻¹, so their waves should have more magnetic than acoustic energy. By comparing Figs. 1 and 2 from Soler et al. (2013b), we see that the wave dissipation is larger for β_plasma ≪ 1 regions than for β_plasma ≫ 1 regions. In addition, the works conducted by Soler et al. (2015a), Popescu Braileanu & Keppens (2021) or Cally & Gómez-Míguez (2023) have shown that if a wave excites efficiently the magnetic field, such as in the case of fast and slow waves when they are magnetic (that is, when β ≪ 1 and β ≫ 1, respectively), or in the case of Alfvén waves, then charge-neutral collisions efficiently damps the wave.

After we analysed the wave energy for each model, we compare the total wave energy from each one. In Fig. 9, we show the difference in wave flux between the models. We observe that the envelope and peak position of the flux differences reminds the behaviour of the heating terms (see Fig. 4, or Fig. 6). We highlight this fact here since it points to a physical origin for the total flux differences. Another important aspect is that the flux differences are mostly positive along the atmosphere, namely, the wave in the 1F model is carrying more energy to the higher layers than in the 2F model. The 2F model describes the dynamics of two fluids, whereas the 1F model restricts itself to the centre of mass of the system, so we would expect the energy of the 2F model to be equal or overcome the 1F energy. Therefore, the result can be explained by considering that the 2F model is more efficient at dissipating energy at deep layers than the 1F model. This view also aligns with the temperature increase differences shown in Fig. 4, where we have noted that the 2F model increases the temperature more than the 1F model before the heating peak.

Fig. 9. Residuals of the total wave energy flux (1F–2F). The results are computed for the different wave periods considered in this work. The vertical dashed line indicates the equipartition layer βplasma = 1.

The analysis of the simulations brings us to some key points about wave heating through charge-neutral collisions. Commonly in literature (i.e. Zaqarashvili et al. 2011a; Soler et al. 2013a,b; Khomenko et al. 2014) it is mentioned the ratio, ω/ν, with ν the collision frequency of the considered charge-neutral model, as the criterion for considering the partially ionised effects. However, in line with the findings of Popescu Braileanu et al. (2019b) or Zhang et al. (2021), as well as this work, the most significant temperature increases are not driven by the short period waves. In those set-ups, the waves are launched at the bottom of the atmosphere and propagate vertically. Then, the long-period waves can transport more energy than short-period waves to the high layers, where the coupling between the species weakens. Therefore, the ratio ω/ν is a good reference for a local analysis, but also the capacity of the wave for transporting energy should be taken into account. We also find that when the magnetic energy of the wave increases, the charge-neutral dissipation becomes more efficient, in agreement with the studies by Forteza et al. (2006), Soler et al. (2015a), Popescu Braileanu & Keppens (2021) or Cally & Gómez-Míguez (2023). Finally, we propose the efficiency of 2F model in dissipating energy at deep layers as the explanation for the discrepancy on the temperature increase. However, the physical assumptions that lead to a different dissipation between models are still not clear and will drive the discussion of the following sections.

4. Vertical propagation of MA waves in a stratified atmosphere with a horizontal and homogeneous magnetic field. Analytic approach

In light of the numerical experiments, we went on to apply an analytical approach to deepen into the solutions of each model and to better understand the discrepancies between the models. As we applied the same methodology for both models, we present below all the partial results in parallel. Our first step involved deriving the linear wave equations. We considered the case of vertical propagation (z-axis) in a stratified atmosphere, with a horizontal magnetic field in the x-axis. As usual, we split the variables, say f, into equilibrium and perturbation parts, f = f₀ + f₁, such as |f₁|/f₀ ≪ 1, where the subindex ‘0’ stands for the equilibrium and the subindex ‘1’ stands for the perturbation. We considered a static equilibrium, so u_α, z, 0 = 0 and |u_α, z, 1| much smaller than the phase speed. Subsequently, we removed the subindex ‘1’ from the velocity perturbation. Keeping only the terms up to the first order, the set of 2F linear equations is (see e.g. Popescu Braileanu et al. 2019b):

$$\begin{array}{cc}& \frac{\partial {\rho }_{\mathrm{n},1}}{\partial t}+\frac{\partial }{\partial z}({\rho }_{\mathrm{n},0}{u}_{\mathrm{n},\mathrm{z}})=0,\hfill \end{array}$$(37a)

$$\begin{array}{cc}& \frac{\partial {\rho }_{\mathrm{c},1}}{\partial t}+\frac{\partial }{\partial z}({\rho }_{\mathrm{c},0}{u}_{\mathrm{c},\mathrm{z}})=0,\hfill \end{array}$$(37b)

$$\begin{array}{cc}& {\rho }_{n,0}\frac{\partial u_{\mathrm{n},\mathrm{z}}}{\partial t}+\frac{\partial p_{\mathrm{n},1}}{\partial z}=-{\rho }_{\mathrm{n},1}g+{\alpha }_0{\rho }_{\mathrm{c},0}{\rho }_{\mathrm{n},0}(u_{\mathrm{c},\mathrm{z}}-u_{\mathrm{n},\mathrm{z}}),\hfill \end{array}$$(37c)

$$\begin{array}{cc}& {\rho }_{\mathrm{c},0}\frac{\partial u_{\mathrm{c},\mathrm{z}}}{\partial t}+\frac{\partial p_{\mathrm{c},1}}{\partial z}=-{\rho }_{\mathrm{c},1}g-{\alpha }_0{\rho }_{\mathrm{c},0}{\rho }_{\mathrm{n},0}(u_{\mathrm{c},\mathrm{z}}-u_{\mathrm{n},\mathrm{z}})-\frac{B_{\mathrm{x},0}}{{\mu }_0}\frac{\partial B_{\mathrm{x},1}}{\partial z},\hfill \end{array}$$(37d)

$$\begin{array}{cc}& \frac{\partial p_{\mathrm{n},1}}{\partial t}+u_{\mathrm{n},\mathrm{z}}\frac{d{p}_{\mathrm{n},0}}{\mathit{dz}}-c_{\mathrm{sn}}^2(\frac{\partial {\rho }_{\mathrm{n},1}}{\partial t}+u_{\mathrm{n},\mathrm{z}}\frac{d{\rho }_{\mathrm{n},0}}{\mathit{dz}})=0,\hfill \end{array}$$(37e)

$$\begin{array}{cc}& \frac{\partial p_{\mathrm{c},1}}{\partial t}+u_{\mathrm{c},\mathrm{z}}\frac{d{p}_{\mathrm{c},0}}{\mathit{dz}}-c_{\mathrm{sc}}^2(\frac{\partial {\rho }_{\mathrm{c},1}}{\partial t}+u_{\mathrm{c},\mathrm{z}}\frac{d{\rho }_{\mathrm{c},0}}{\mathit{dz}})=0,\hfill \end{array}$$(37f)

$$\begin{array}{cc}& \frac{\partial B_{\mathrm{x},1}}{\partial t}=-B_{\mathrm{x},0}\frac{\partial u_{\mathrm{c},\mathrm{z}}}{\partial z}·\hfill \end{array}$$(37g)

Similarly, the 1F set of linear equations is

$$\begin{array}{cc}& \frac{\partial {\rho }_1}{\partial t}+\frac{\partial }{\partial z}({\rho }_0{u}_{\mathrm{z}})=0,\hfill \end{array}$$(38a)

$$\begin{array}{cc}& {\rho }_0\frac{\partial u_{\mathrm{z}}}{\partial t}+\frac{\partial p_1}{\partial z}=-{\rho }_{1}g-\frac{B_{\mathrm{x},0}}{{\mu }_0}\frac{\partial B_{\mathrm{x},1}}{\partial z},\hfill \end{array}$$(38b)

$$\begin{array}{cc}& \frac{\partial p_1}{\partial t}+u_{\mathrm{z}}\frac{d{p}_0}{\mathit{dz}}-c_s^2(\frac{\partial {\rho }_1}{\partial t}+u_z\frac{d{\rho }_0}{\mathit{dz}})=0,\hfill \end{array}$$(38c)

$$\begin{array}{cc}& \frac{\partial B_{\mathrm{x},1}}{\partial t}=-B_{\mathrm{x},0}\frac{\partial u_{\mathrm{z}}}{\partial z}+c_{\mathrm{A}}^2{\rho }_0\frac{\partial }{\partial z}\left({\stackrel{\sim }{\eta }}_{\mathrm{A},0}\frac{\partial B_{\mathrm{x},1}}{\partial z}\right),\hfill \end{array}$$(38d)

where the only non-ideal contribution comes from the ambipolar diffusion in the induction equation.

By combining the linear equations, we obtain a set of two second order linear equations that involves the velocities of each fluid. For the 2F model, we find

$$\begin{array}{cc}\hfill \frac{{\partial }^2{u}_{\mathrm{n},\mathrm{z}}}{\partial t^2}& =c_{\mathrm{sn}}^2\frac{{\partial }^2{u}_{\mathrm{n},\mathrm{z}}}{\partial z^2}+\frac{1}{{\rho }_{\mathrm{n},0}}\frac{d}{\mathit{dz}}(c_{\mathrm{sn}}^2{\rho }_{n,0})\frac{\partial u_{\mathrm{n},\mathrm{z}}}{\partial z}\hfill \\ \hfill & +{\alpha }_0{\rho }_{\mathrm{c},0}\frac{\partial }{\partial t}(u_{\mathrm{c},\mathrm{z}}-u_{\mathrm{n},\mathrm{z}}),\hfill \end{array}$$(39a)

$$\begin{array}{cc}\hfill \frac{{\partial }^2{u}_{\mathrm{c},\mathrm{z}}}{\partial t^2}& =(c_{\mathrm{sc}}^2+c_{\mathrm{Ac}}^2)\frac{{\partial }^2{u}_{\mathrm{c},\mathrm{z}}}{\partial z^2}+\frac{1}{{\rho }_{\mathrm{c},0}}\frac{d}{\mathit{dz}}[(c_{\mathrm{sc}}^2+c_{\mathrm{Ac}}^2){\rho }_{\mathrm{c},0}]\frac{\partial u_{\mathrm{c},\mathrm{z}}}{\partial z}\hfill \\ \hfill & -{\alpha }_0{\rho }_{\mathrm{n},0}\frac{\partial }{\partial t}(u_{\mathrm{c},\mathrm{z}}-u_{\mathrm{n},\mathrm{z}}),\hfill \end{array}$$(39b)

and, for the 1F model

$$\begin{array}{cc}\hfill \frac{{\partial }^2{u}_{\mathrm{z}}}{\partial t^2}& =c_{\mathrm{s}}^2\frac{{\partial }^2{u}_{\mathrm{z}}}{\partial z^2}+\frac{1}{{\rho }_0}\frac{d}{\mathit{dz}}(c_{\mathrm{s}}^2{\rho }_0)\frac{\partial u_{\mathrm{z}}}{\partial z}+\frac{1}{{\rho }_0{\stackrel{\sim }{\eta }}_{\mathrm{A},0}}\frac{\partial }{\partial t}(u_{\mathrm{Amb},\mathrm{z}}-u_{\mathrm{z}}),\hfill \end{array}$$(40a)

$$\begin{array}{cc}\hfill 0& =c_{\mathrm{A}}^2\frac{{\partial }^2{u}_{\mathrm{Amb},\mathrm{z}}}{\partial z^2}-\frac{1}{{\rho }_0{\stackrel{\sim }{\eta }}_{\mathrm{A},0}}\frac{\partial }{\partial t}(u_{\mathrm{Amb},\mathrm{z}}-u_{\mathrm{z}}),\hfill \end{array}$$(40b)

where we introduce the charge velocity due to ambipolar diffusion, u_Amb, z, as

$$\begin{array}{c}\hfill u_{\mathrm{Amb},\mathrm{z}}=u_{\mathrm{z}}-{\stackrel{\sim }{\eta }}_{\mathrm{A},0}\frac{B_{\mathrm{x},0}}{{\mu }_0}\frac{\partial B_{\mathrm{x},1}}{\partial z}·\end{array}$$(41)

This definition follows from introducing w_1F (Eq. (23)) in relation (18), only accounting for the magnetic term. In the 2F model, we find wave equations for the velocities of neutrals and charges. Conversely, a wave equation couples with a diffusion equation in the 1F model. We also note that a term proportional to $\frac{\partial u_{\mathrm{Amb},\mathrm{z}}}{\partial z}$, which accounts for the stratification, is absent from the equation of u_Amb, z because B_x, 0 is uniform.

Across this work, we have denoted the decay of the wave amplitude as damping or dissipation. In absence of other effects, wave dissipation implies that the wave amplitude decreases strictly because the wave is releasing energy to the environment. Note however that there is no dissipation term in the linear energy equations as that term is of the second order, therefore, heating cannot be properly captured on the linear approach (see e.g. Cally & Gómez-Míguez 2023; Cally 2023). In the case of the stratified atmosphere, velocity waves increase their amplitude but not their energy flux due to the gravitational stratification. In terms of velocity amplitude, the growth by stratification compensates for the effects of the dissipation when dissipation is still weak. Consequently, Figs. 7 and 8 show that wave energy flux is constant at deep layers, and above them the energy flux drops down. The velocity waves we are studying propagate to layers with decreasing density, increasing their amplitude, and to layers of decreasing collision frequency, enhancing the decoupling effects. The latter stands for a propagating wave such ω < αρ₀. Propagating magneto-acoustic waves exhibits a maximum dissipation for ω ≈ αρ₀ and it decreases for larger frequencies (see discussions by Soler et al. 2013a; Popescu Braileanu et al. 2019b). For the waves we study, both effects are entangled, so it is not straightforward to talk about dissipation properly. Equations (39a), (39b), (40a) and (40b) can be solved by variable separation. If the background variables are constant, then the wave variables admit a solution of the type $f(z,t)=\stackrel{\sim }{f}(z)e^{-i\omega t}$. The term proportional to the first spatial derivative is the one that drives the change of the wave amplitude along the atmosphere. To decouple this effect from the dissipation, we introduce the wave variables (see similar derivations by Lamb 1932; Deubner & Gough 1984):

$$\begin{array}{cc}\hfill {\psi }_{\mathrm{n},\mathrm{z}}& ={({\rho }_{\mathrm{n},0}{c}_{\mathrm{sn}}^2)}^{\frac{1}{2}}{\stackrel{\sim }{u}}_{\mathrm{n},\mathrm{z}}=a_{\mathrm{n}}{\stackrel{\sim }{u}}_{\mathrm{n},\mathrm{z}},\hfill \end{array}$$(42)

$$\begin{array}{cc}\hfill {\psi }_{\mathrm{c},\mathrm{z}}& ={({\rho }_{\mathrm{c},0}(c_{\mathrm{sc}}^2+c_{\mathrm{A},\mathrm{c}}^2))}^{\frac{1}{2}}{\stackrel{\sim }{u}}_{\mathrm{c},\mathrm{z}}=a_{\mathrm{c}}{\stackrel{\sim }{u}}_{\mathrm{c},\mathrm{z}},\hfill \end{array}$$(43)

$$\begin{array}{cc}\hfill {\psi }_{\mathrm{z}}& ={({\rho }_0{c}_{\mathrm{s}}^2)}^{\frac{1}{2}}{\stackrel{\sim }{u}}_{\mathrm{z}}=a_{\mathrm{u}}{\stackrel{\sim }{u}}_{\mathrm{z}},\hfill \end{array}$$(44)

$$\begin{array}{cc}\hfill {\psi }_{\mathrm{Amb},\mathrm{z}}& ={({\rho }_0{c}_{\mathrm{A}}^2)}^{\frac{1}{2}}{\stackrel{\sim }{u}}_{\mathrm{Amb},\mathrm{z}}=a_{\mathrm{Amb}}{\stackrel{\sim }{u}}_{\mathrm{Amb},\mathrm{z}}.\hfill \end{array}$$(45)

The transformations are introduced in Eqs. (39a), (39b), (40a), and (40b), taking into account that the proportional factors are not a function of the time. Now, the ψ functions only depend on z and the time derivatives can be removed from the wave equations accordingly to the sinusoidal solution. Therefore, the equations for the 2F spatial variables ψ read as (in 2F),

$$\begin{array}{cc}\hfill \frac{d^2{\psi }_{\mathrm{n},\mathrm{z}}}{d{z}^2}+\frac{{\omega }^2}{c_{\mathrm{sn}}^2}{\psi }_{\mathrm{n},\mathrm{z}}& =i\omega {\alpha }_0\frac{{\rho }_{\mathrm{c},0}}{c_{\mathrm{sn}}^2}(\frac{a_{\mathrm{n}}}{a_{\mathrm{c}}}{\psi }_{\mathrm{c},\mathrm{z}}-{\psi }_{\mathrm{n},\mathrm{z}}),\hfill \end{array}$$(46a)

$$\begin{array}{cc}\hfill \frac{d^2{\psi }_{\mathrm{c},\mathrm{z}}}{d{z}^2}+\frac{{\omega }^2}{c_{\mathrm{sc}}^2+c_{\mathrm{Ac}}^2}{\psi }_{\mathrm{c},\mathrm{z}}& =-i\omega {\alpha }_0\frac{{\rho }_{n,0}}{c_{\mathrm{sc}}^2+c_{\mathrm{Ac}}^2}({\psi }_{\mathrm{c},\mathrm{z}}-\frac{a_{\mathrm{c}}}{a_{\mathrm{n}}}{\psi }_{\mathrm{n},\mathrm{z}}),\hfill \end{array}$$(46b)

and those for the 1F model (see Eqs. (40a) and (40b))

$$\begin{array}{cc}\hfill \frac{d^2{\psi }_{\mathrm{u},\mathrm{z}}}{d{z}^2}+\frac{{\omega }^2}{c_{\mathrm{s}}^2}{\psi }_{\mathrm{u},\mathrm{z}}& =\frac{i\omega }{{\rho }_0{c}_{\mathrm{s}}^2{\eta }_{\mathrm{A},0}}(\frac{a_{\mathrm{u}}}{a_{\mathrm{Amb}}}{\psi }_{\mathrm{Amb},\mathrm{z}}-{\psi }_{\mathrm{u},\mathrm{z}}),\hfill \end{array}$$(47a)

$$\begin{array}{cc}\hfill \frac{d^2{\psi }_{\mathrm{Amb},\mathrm{z}}}{d{z}^2}& =-\frac{i\omega }{c_{\mathrm{A}}^2{\rho }_0{\eta }_{\mathrm{A},0}}({\psi }_{\mathrm{Amb},\mathrm{z}}-\frac{a_{\mathrm{Amb}}}{a_{\mathrm{u}}}{\psi }_{\mathrm{u},\mathrm{z}}).\hfill \end{array}$$(47b)

We neglected the terms that involve the derivatives of a_α terms, with α ∈ {n, c, u, Amb}, consistently with the condition that the wavelength is much smaller than the pressure scale height.

We used the Eikonal approach to approximate the wave equation solution, assuming that the wave amplitude changes on a scale larger than its wavelength (see e.g. Weinberg 1962; de Pontieu & Haerendel 1998; Tracy et al. 2014). This means that a solution of the form ψ = A(z)e^iθe(z) can be proposed, where ψ is the vector for the wave solutions (e.g. ψ = (ψ_n, ψ_c) in the 2F model), A(z) is the wave amplitude, θ = ∫_z0^zk_z(z′)dz′ is the spatial wave phase, with k_z the vertical wave number, and e is the polarisation vector. By neglecting spatial derivatives on k_z, A and e, (see e.g. de Pontieu & Haerendel 1998), the set of differential equations reduces to a set of polynomial equations of the form

$$\begin{array}{c}\hfill {\widehat{D}}_{2\mathrm{F}}·{\mathbf{e}}_{2\mathrm{F}}=0;{\widehat{D}}_{1\mathrm{F}}·{\mathbf{e}}_{1\mathrm{F}}=0,\end{array}$$(48)

where

$$\begin{array}{cc}\hfill {\widehat{D}}_{2F}& =\left(\begin{array}{cc}{k}_{\mathrm{z}}^2-(\frac{{\omega }^2}{c_{\mathrm{sn}}^2}+i\omega {\alpha }_0\frac{{\rho }_{\mathrm{c},0}}{c_{\mathrm{sn}}^2})& i\omega {\alpha }_0\frac{{\rho }_{\mathrm{c},0}}{c_{\mathrm{sn}}^2}\frac{a_{\mathrm{n}}}{a_{\mathrm{c}}}\\ i\omega {\alpha }_0\frac{{\rho }_{\mathrm{n},0}}{c_{\mathrm{sc}}^2+c_{\mathrm{Ac}}^2}\frac{a_{\mathrm{c}}}{a_{\mathrm{n}}}& k_{\mathrm{z}}^2-(\frac{{\omega }^2}{c_{\mathrm{sc}}^2+c_{\mathrm{Ac}}^2}+i\omega {\alpha }_0\frac{{\rho }_{\mathrm{n},0}}{c_{\mathrm{sc}}^2+c_{\mathrm{Ac}}^2})\end{array}\right),\hfill \end{array}$$(49)

$$\begin{array}{cc}\hfill {\widehat{D}}_{1F}& =\left(\begin{array}{cc}{k}_{\mathrm{z}}^2-(\frac{{\omega }^2}{c_{\mathrm{s}}^2}+i\frac{\omega }{c_{\mathrm{s}}^2{\rho }_0{\stackrel{\sim }{\eta }}_{\mathrm{A},0}})& i\frac{\omega }{c_{\mathrm{s}}^2{\rho }_0{\eta }_{\mathrm{A},0}}\frac{a_{\mathrm{u}}}{a_{\mathrm{Amb}}}\\ i\frac{\omega }{c_{\mathrm{A}}^2{\rho }_0{\stackrel{\sim }{\eta }}_{\mathrm{A},0}}\frac{a_{\mathrm{Amb}}}{a_{\mathrm{u}}}& k_{\mathrm{z}}^2-i\frac{\omega }{c_{\mathrm{A}}^2{\rho }_0{\eta }_{\mathrm{A},0}}\end{array}\right),\hfill \end{array}$$(50)

are the dispersion matrices and

$$\begin{array}{c}\hfill \mathbf{e}=\frac{1}{\sqrt{{\left|\frac{D_{12}}{D_{11}}\right|}^2+1}}(-\frac{D_{12}}{D_{11}},1)\end{array}$$(51)

is the polarisation vector, which is defined as the unit null vector of the dispersion matrix.

The determinant of the dispersion matrix must be equal to zero to obtain the non-trivial solution of the system of equations. This condition leads to the well-known dispersion relations (see (Kumar & Roberts 2003; Zaqarashvili et al. 2011a; Soler et al. 2013b; Popescu Braileanu et al. 2019b) for the 2F case and (Chapman & Cowling 1970; Forteza et al. 2006; Ballester et al. 2018) for the 1F one), which are, in each case,

$$\begin{array}{cc}\hfill det({\widehat{D}}_{2\mathrm{F}})& =k_{\mathrm{z}}^4-[\frac{{\omega }^2}{c_{\mathrm{sn}}^2}+\frac{{\omega }^2}{c_{\mathrm{sc}}^2+c_{\mathrm{Ac}}^2}+i\omega {\alpha }_0(\frac{{\rho }_{\mathrm{c},0}}{c_{\mathrm{sn}}^2}+\frac{{\rho }_{\mathrm{n},0}}{c_{\mathrm{sc}}^2+c_{\mathrm{Ac}}^2})]k_{\mathrm{z}}^2\hfill \\ \hfill & +\frac{{\omega }^4}{c_{\mathrm{sn}}^2(c_{\mathrm{sc}}^2+c_{\mathrm{Ac}}^2)}+i{\omega }^3{\alpha }_0\left(\frac{{\rho }_{\mathrm{n},0}+{\rho }_{\mathrm{c},0}}{c_{\mathrm{sn}}^2(c_{\mathrm{sc}}^2+c_{\mathrm{Ac}}^2)}\right)=0,\hfill \end{array}$$(52)

$$\begin{array}{cc}\hfill det({\widehat{D}}_{1\mathrm{F}})& =k_{\mathrm{z}}^4-[\frac{{\omega }^2}{c_{\mathrm{s}}^2}+i\frac{\omega }{{\rho }_0{\stackrel{\sim }{\eta }}_{\mathrm{A},0}}(\frac{1}{c_{\mathrm{s}}^2}+\frac{1}{c_{\mathrm{A}}^2})]k_{\mathrm{z}}^2+i\frac{{\omega }^3}{c_{\mathrm{s}}^2{c}_{\mathrm{A}}^2{\rho }_0{\stackrel{\sim }{\eta }}_{\mathrm{A},0}}=0.\hfill \end{array}$$(53)

Both equations are 4th order in wave number (k_z), leading to two different pairs of normal modes, each pair containing an upward and a downward propagating mode. We observe that $det({\widehat{D}}_{2\mathrm{F}})$ contains a term proportional to ω⁴ that is missing in the $det({\widehat{D}}_{1\mathrm{F}})$ expression. From Eq. (52), the ratio of the ω⁴-proportional term to the ω³-proportional one, in absolute value, is ω/(ρ₀α₀). This ratio is nothing but the estimator for the contribution of the inertial terms in the drift velocity equation (see Eq. (22)). Therefore, the absence of the ω⁴ term is the consequence of neglecting inertial terms in Eq. (23), as noted by Zaqarashvili et al. (2011a) for the case of Alfvén waves.

Working in terms of the wave functions described by Eqs. (42), (43), (44), and (45), rather than the model variables (u_n, u_c, u, B_x) offers the advantage of symmetry between upward and downward modes, and the imaginary part of the wave number is directly linked to the collisional dissipation. In Fig. 10, we show the solutions of the dispersion relations given by Eqs. (52) and (53) for our model atmosphere.

Fig. 10. Vertical wave numbers for the 2F model as functions of height for all wave periods studied in this work. From top to bottom, the wave periods are 1 s, 3 s, and 5 s. Blue and green solid lines: upward and downward fast mode. Red and black solid lines: Upward and downward normal mode (different from the fast mode). The left column represents the real part and the right represents the imaginary part. The fast mode lines are multiplied by an arbitrary factor for the representation. The vertical dashed line indicates the equipartition layer βplasma = 1.

One pair of solutions has a similar behaviour as the classic fast magneto-acoustic mode, but it dissipates energy because of charge-neutral collisions (mode 4 travels upwards and mode 3, downwards). As we have seen before, at deep layers the fast mode behaves as an acoustic wave, travelling at the speed c_s ∼ c_sn, since the atmosphere is weakly ionised. Once the wave crosses the equipartition layer, it behaves as a magnetic wave, travelling with a speed that tends to c_A when reaching layers where β_plasma ≪ 1. The other mode (mode 2 traveling upwards and mode 1, downwards) is highly dissipated at the deep layers as the real and imaginary parts of k_z have similar values. These solutions resemble those presented in Fig. 13 of Popescu Braileanu et al. (2019b), but here we have included all the wave periods explored in this work.

For model comparison, we show in Fig. 11 (black solid line) the relative difference between the 1F and 2F wave numbers with respect to the value computed for the 2F model, for both the real and imaginary parts. At this point of the discussion, we focus on the black curve, letting the others for Sect. 5.1. We only excite the fast mode in the numerical experiments, so we focus the discussion on this mode. Only the upwards mode is represented. Differences in the real part of k_z are negligible at deep layers (of the order of 0.1%), reaching a maximum of 1% in the higher layers for the 1 s wave period. As expected, those differences increase with the wave frequency. However, the imaginary part of k_z shows more significant differences at deep and high layers. We note that the differences at the high layers increase with the wave frequency, but those at deep layers seem to be unaffected by the wave period. This result shows that the waves are more damped in the 2F model at deep layers, but the opposite occurs in the high layers. This picture matches qualitatively with the temperature increase residuals shown in Fig. 4. Even though Fig. 6 informs us that the non-linear effects play an important role, Fig. 11 gives us a chance to figure out an explanation for model discrepancies using the linear theory.

Fig. 11. Difference between the 2F and the 1F wave numbers relative to the 2F wave number. Left column: Calculations for the real part of kz; right column: same calculations, but for the imaginary part of kz. Only the upwards fast mode solution is represented. Solid black line: differences computed asumming the 1F model described in Section 2.2; dashed red line: differences computed extending the 1F for accounting the contributions of qw and G; dashed-dotted red line: differences computed extending the 1F for accounting the contributions of qw alone and G alone. The vertical dashed line indicates the equipartition layer βplasma = 1.

A useful parameter to put into context the importance of the differences in k_z that we have obtained is the quality factor of the waves, commonly defined as

$$\begin{array}{c}\hfill Q_{k_{\mathrm{z}}}=\frac{1}{2}\left|\frac{\mathfrak{R}(k_{\mathrm{z}})}{\mathfrak{I}(k_{\mathrm{z}})}\right|.\end{array}$$(54)

If $Q_{k_{\mathrm{z}}}=\frac{1}{2}$, the wave is effectively damped in a distance of a wavelength. For $Q_{k_{\mathrm{z}}}>\frac{1}{2}$, the wave dissipates energy in scales larger than its wavelength and for $Q_{k_{\mathrm{z}}}<\frac{1}{2}$, it is efficiently dissipated before reaching the scale of its wavelength. Fig. 12 shows that below 1000 km, dissipating the wave energy requires more than 100 periods, indicating a weak dissipation. However, non-linear effects become important for the 1 s and 3 s period waves at that height, enhancing energy dissipation, as shown in Figs. 6, 7 and 8. As the waves reach layers of β_plasma < 1, the wavelength is lengthened and the imaginary part of k_z increases, so we get Q ≈ 1 at the high layers (see Fig. 10). From Fig. 11, we see that the more significant differences appear at the bottom of the atmosphere, where the waves propagate almost without dissipation or at the top layers, where the waves arrive with low energy. Therefore, the differences in wave flux shown in Fig. 9 may be caused by deviations at the intermediate layers.

Fig. 12. Quality factors for the simulated wave periods for the 2F model. The horizontal grey line indicates the value $Q_{k_{\mathrm{z}}}^{2\mathrm{F}}=\frac{1}{2}$. The vertical dashed line indicates the equipartition layer βplasma = 1.

5. Model discrepancies

In the previous sections, we compare a purely hydrogen 2F model with a standard 1F model with ambipolar diffusion, whose parameters ${\xi }_{\mathrm{n}},\stackrel{\sim }{\mu },{\stackrel{\sim }{\eta }}_{\mathrm{A}}$ are introduced consistently for such plasma. After examining numerical simulations and analytical solutions, we found deviations in the imaginary part of k_z between models to explain why the wave in the 1F model is upwardly transporting more energy than in the 2F model. Next, we tried to determine the origin of those differences by looking into the model construction. As the main assumptions of the 1F model are listed in Section 2.2, here we delve into the variable definitions.

While the frame of reference in NF models refers, say, to the laboratory frame of reference, the 1F variables refer to the centre of mass frame of reference, so their definitions depend on the relative velocities of each component to the centre of mass velocity, w_α. In Appendix B, we list the definitions of the 1F variables in terms of the 2F variables. For further details, we refer to the works by Braginskii (1965), Schunk (1977) or Khomenko et al. (2014), among others. The inspection of the definitions bring out that most of them involve terms proportional to ρ_αw_α² or ρ_αw_α³, which can be neglected in our study (see discussion above Eq. (34)). However, the 1F heat flux (see Eq. (B.1d)) includes two terms proportional to w_α and to the thermal pressure. The latter was also noted by Schunk (1977), who suggested checking these terms when comparing models with different variable definitions in the velocity frame of reference. The 1F heat flux components proportional to w_α only appear when it has been derived from a NF model closed without including a heat flux for each species, as the 2F model discussed here. In such situation, as shown in Appendix C, the terms regarding the 1F heat flux definition appear, even though there is not any heat flux present in the NF model. From Eqs. (B.1d), (17a) and (17b), we find this 1F heat flux to be

$$\begin{array}{c}\hfill {\mathbf{q}}_{\mathrm{w}}\frac{\gamma }{\gamma -1}{\displaystyle \sum _{\{\mathrm{n},\mathrm{c}\}}}{p}_{\alpha }{\mathbf{w}}_{\alpha }\approx \frac{\gamma }{\gamma -1}{\xi }_{\mathrm{n}}(1-{\xi }_{\mathrm{n}})\stackrel{\sim }{\mu }p\mathbf{w},\end{array}$$(55)

in absence of viscosity. We point to the second expression when referring to q_w, which is derived from the previous one under the assumption of T_n = T_c, as in the case of G (see Eq. (25)). In addition, as for the G term, it vanishes for the high and low ionisation limits. Both S^* (or Q) and q_w are of the order 𝒪(w), so the last one seems a valuable candidate to take into account for explaining the model discrepancies.

5.1. A look into the neglected contributions

Following our discussion of the main approximations considered in the 1F modelling, we examine the ones we consider that mainly influence the discrepancies between models: the contributions of G, q_w and the inertial terms. For the first two cases, Eqs. (25) and (55) allow us to incorporate these effects in the 1F model described in Section 2.2. For the case of the inertial terms, it is not straightforward to account for them easily in the 1F model, but we can use Eq. (22) for estimating them.

We first considered the contributions of G and q_w to the dispersion relation. For the sake of simplicity, we directly neglected the terms involving the derivatives of the background quantities, which is reasonable for the range of frequencies of this study. The mass and momentum equations remain as before. We leave the details about how to extend the energy and induction equation to Appendix C. From the set of Eqs. (C.3a), (C.3b), (C.3c), and (C.3d), we obtain the dispersion relation

$$\begin{array}{c}\hfill k_{\mathrm{z}}^4-{(1+{\xi }_{\mathrm{c},0})}^2[\frac{{\omega }^2}{c_{\mathrm{s}}^2}+\frac{{\omega }^2{\xi }_{\mathrm{c},0}^2{\stackrel{\sim }{\mu }}_0^2}{c_{\mathrm{A}}^2}+i\frac{\omega }{{\rho }_0{\stackrel{\sim }{\eta }}_{\mathrm{A},0}}(\frac{1}{c_{\mathrm{s}}^2}+\frac{1}{c_{\mathrm{A}}^2})]k_{\mathrm{z}}^2\\ \hfill +i{(1+{\xi }_{\mathrm{c},0})}^2\frac{{\omega }^3}{c_{\mathrm{s}}^2{c}_{\mathrm{A}}^2{\rho }_0{\stackrel{\sim }{\eta }}_{\mathrm{A},0}}=0.\end{array}$$(56)

As expected, G and q_w slightly modify the expression derived in Eq. (53). We solve Eq. (56) for the upwards fast mode and represent it in Fig. 11 in red dashed lines. We refer to this model as the extended 1F model. We find that the new dispersion relation can solve the discrepancies between models at deep and intermediate layers. As we have noted before (see discussion in Section 4 about Fig. 11), those differences look weakly dependent on the wave period. For the model atmosphere, we have $\frac{{\xi }_{\mathrm{c},0}^2}{c_{\mathrm{A}}^2}=\frac{1}{c_{\mathrm{A},\mathrm{c}}^2}\ll \frac{1}{c_{\mathrm{s}}^2}$ and the terms proportional to ω² in the coefficient of k_z² in Eq. (56) reduce to the ones in Eq. (53). Therefore, the modifications introduced by G and q_w (proportional to 1 + ξ_c, 0) are appropriate for explaining those differences. As closure, if one only considers the term G or the term q_w (only the ambipolar contribution) alone, it follows

$$\begin{array}{c}\hfill k_{\mathrm{z}}^4-(1+{\xi }_{\mathrm{c},0})[\frac{{\omega }^2}{c_{\mathrm{s}}^2}+i\frac{\omega }{{\rho }_0{\stackrel{\sim }{\eta }}_{\mathrm{A},0}}(\frac{1}{c_{\mathrm{s}}^2}+\frac{1}{c_{\mathrm{A}}^2})]k_{\mathrm{z}}^2\\ \hfill +i(1+{\xi }_{\mathrm{c},0})\frac{{\omega }^3}{c_{\mathrm{s}}^2{c}_{\mathrm{A}}^2{\rho }_0{\stackrel{\sim }{\eta }}_{\mathrm{A},0}}=0.\end{array}$$(57)

This dispersion relation is a particular case of the one derived by Forteza et al. (2006), when setting k_x = 0 and in absence of collisions with electrons. We also represent this case in Fig. 11, finding that both terms are needed for the correct match between the imaginary part of the wave numbers.

After adding the corrections of our extended 1F model, the period dependent discrepancies at the high layers remain. In Fig. 13 we represent the estimator for the inertial terms for all wave periods computed after Eq. (22). We find them to be negligible at the deep layers, where the collision frequency is large, but they become more relevant at the high layers. By comparing ϖ_inertia with the curves for the extended 1F model in Fig. 11, we find the same trend with height and period, and the values are comparable. As mentioned, we cannot provide a direct proof confirming that inertial terms produce those differences at high layers, however, Fig. 11 tentatively points to such conclusion.

Fig. 13. Estimator of the impact of inertial terms for the three wave periods considered in this work as function of the height. The vertical dashed line indicates the equipartition layer βplasma = 1.

It must be noted that the above discussion is limited to the model atmosphere considered in the present work, which has a rather low ionisation fraction. For high ionisation fractions, it can be expected that both magnetic and thermal pressure decoupling mechanisms vanish (see Eq. (20) or Eq. (29) in Forteza et al. 2006) and the heating is significantly reduced. Therefore, the discrepancies between 2F and 1F models would be also reduced in absolute terms. In the intermediate situation, when neutrals and charges have similar numbers, both decoupling mechanisms become of the same order of magnitude (see discussion by Pinto et al. 2008). The drift velocities are expected to be smaller than in the weakly ionised atmosphere, but, nevertheless, non negligible.

5.2. Implications on the non-linear regime

We present an argument for solving the discrepancies between models in the linear regime above. Next, we take a look back to the simulations and the non-linear regime.

We focus again on the internal energy equations. We established that the main discrepancies between models are driven by the intermediate layers, so then we considered the approximation Eq. (23) to the drift velocity and leave the inertial terms aside. Then, we took Eq. (28) for the 2F model and Q = J ⋅ E^* for the 1F model, so that

$$\begin{array}{cc}\hfill Q_{\mathrm{n}}+Q_{\mathrm{c}}& =\frac{{\stackrel{\sim }{\eta }}_{\mathrm{A}}}{{\xi }_{\mathrm{n}}^2}[{\xi }_{\mathrm{n}}^2{(\mathrm{\nabla }·{\widehat{p}}_{\mathrm{m}})}^2+2{\xi }_{\mathrm{n}}\mathbf{G}·(\mathrm{\nabla }·{\widehat{p}}_{\mathrm{m}})+G^2],\hfill \end{array}$$(58a)

$$\begin{array}{cc}\hfill Q& =\frac{{\stackrel{\sim }{\eta }}_{\mathrm{A}}}{{\xi }_{\mathrm{n}}^2}[{\xi }_{\mathrm{n}}^2{(\mathrm{\nabla }·{\widehat{p}}_{\mathrm{m}})}^2+{\xi }_{\mathrm{n}}\mathbf{G}·(\mathrm{\nabla }·{\widehat{p}}_{\mathrm{m}})].\hfill \end{array}$$(58b)