© ESO 2019

1. Introduction

Solar chromosphere is a transition layer where the properties of plasma change abruptly from being dominated by gas pressure (as in the photosphere) to being dominated by magnetic field (in corona), and where the collisional coupling of the plasma weakens with height. Several aspects of chromospheric physics make this region one of the most complex for understanding of the solar atmosphere.

It is well known that the chromosphere is not in local thermodynamic equilibrium (LTE). Athay & House (1962) measured and compared intensities in strong emission lines of H, He I, He II, and Ca II, and the deduced ratios of the populations of various levels indicated a significant departure from a Boltzmann distribution of the levels studied. In the high chromosphere the ratio of neutral atoms to ions exceeds the LTE values by about 10⁷–10⁸ in the case of He I and by about 10⁶–10⁷ in the case of hydrogen, and there is evidence that significant departures from LTE extend well into the photospheric layers. The formation of many photospheric lines is well described using LTE when radiation is fully coupled to the local physical conditions. However, in the chromosphere when collisional rates are lower, the radiation begins to decouple from the local physical conditions, invalidating the LTE assumption.

Plasma in the chromosphere is only partially ionized. As the density decreases with height, the collision frequencies between different particle species decrease and the magnetic field begins to dominate the dynamics. The neutral species do not feel the magnetic field directly and this causes a partial decoupling between charges and neutrals when the collisional timescales between ionized and neutral atomic species become equal to or larger than the magnetohydrodynamic (MHD) timescales. Therefore, the classical MHD approach is not the best approximation to treat the plasma processes in the chromosphere. While realistic photospheric models based on the MHD approximation provide results that are practically indistinguishable from observations even in strongly magnetized regions (e.g., Rempel et al. 2009), this is not the case for the chromosphere (see e.g., the discussion in Leenaarts 2010). This is because observations do not reach scales where ion-neutral effects can be detected directly in plasma dynamics (see e.g., Fig. 1 in Khomenko et al. 2014) in the photosphere. However, in the chromosphere, these scales become larger and may be reached in some dedicated observations. Indirect ion-neutral decoupling effects such as heating are also more pronounced in the chromosphere (e.g., see Ni et al. 2018).

A suitable alternative to the MHD approach is a two-fluid plasma–neutrals model, numerically implemented in this work. This kind of model is most appropriate for high-density plasmas such as in the chromosphere, where kinetic models cannot be efficiently applied. Kinetic models have been successfully applied to the corona and the solar wind, where much lower particle density makes the problem more computationally feasible and the kinetic description more necessary and appropriate. Multi-fluid models preserve a number of important physical effects coming from weak collisional coupling, still allowing for an efficient implementation, as in any fluid-based description. While the motivation for multi-fluid modeling comes mainly from theoretical considerations, in recent years there have been attempts at direct observational detection of multi-fluid effects in the solar plasma through measurements of differences in ion and neutral velocities caused by loss of collisional coupling. By measuring the Doppler shift in ion Fe II and neutral Fe I lines simultaneously over the same volume of plasma, differences between ion and neutral velocities of the Evershed flow have been found by Khomenko et al. (2015) in sunspot penumbra as deep as the photosphere. Similarly, Khomenko et al. (2016) showed that non-negligible differences in He I and Ca II velocities exist in solar prominences when observed with sufficiently high cadence. While the measured neutral and ion velocities were the same in most of the spatial locations most of the time, they were observed to decouple in the presence of strong spatial and temporal gradients, such as wave fronts. However, Anan et al. (2017) concluded that the observed small differences in velocities of H I 397 nm, H I 434 nm, Ca II 397 nm, and Ca II 854 were indications of motion of different components in the prominence along the line of sight, rather than indications of decoupling. Direct detection of the decoupling effects would require the highest spatial and temporal resolution, apart from a careful selection of spectral lines for the analysis. More observational studies in the future will be needed to confirm or discard the detection of the decoupling effects.

Nevertheless, there exist other indirect indications of decoupling such as the those observed by Gilbert et al. (2007). These authors compared He I 1083 nm and Hα data in multiple solar prominences in different phases of their lifecycle, and were able to detect the drainage effects across prominence magnetic field with different timescales for helium and hydrogen atoms. Such drainage was predicted earlier by Gilbert et al. (2002) from simplified semi-analytical calculations. Further indirect evidence of ion-neutral decoupling was reported by de la Cruz Rodríguez & Socas-Navarro (2011) who deduced misalignment in the visible direction of chromospheric fibrils and measured the magnetic field vector. These latter authors noted that most of the fibrils aligned with the magnetic field, although there are a few noteworthy cases where significant misalignments occur, well beyond the observational uncertainty.

Misalignment of chromospheric fibrils and magnetic field was confirmed by Martínez-Sykora et al. (2016) by using advanced radiative MHD simulations, which included the effects of ion-neutral interactions in the single-fluid approximation (via ambipolar diffusion). It was shown that the magnetic field is indeed often not well aligned with chromospheric features. The misalignment occurs where the ambipolar diffusion is large, that is, ion and neutral populations decouple as the ion-neutral collision frequency drops, allowing the field to slip through the neutral population. The conditions for misalignment also require currents perpendicular to the field to be strong, and thermodynamic timescales to be longer than or similar to those of ambipolar diffusion.

The vast majority of studies of ion-neutral effects in solar plasmas address propagation of different types of waves and the development of instabilities in chromospheric and prominence plasma conditions. Zaqarashvili et al. (2011b) and Soler et al. (2013a) performed theoretical studies of the propagation of Alfvén waves in a uniform medium composed of protons and hydrogen atoms which interact by collisions. Refinements of this theory were introduced by Zaqarashvili et al. (2011a) by adding helium atoms, by Zaqarashvili et al. (2013) by considering a nonuniform gravitationally stratified atmosphere, and by Soler et al. (2013b) by applying a flux tube model. It was found that in situations when the neutral-ion collisional frequency is much lower than the wave frequency, the propagation properties of the waves only depend on the properties of the ionized fluid. These properties change with height due to stratification, and the Alfvén wave may become evanescent in some regions of the solar atmosphere. The waves can be damped by neutral-ion collisions, and the existence of neutral helium atoms, alongside with neutral hydrogen, significantly enhances the damping. The damping is most efficient when the wave frequency and the collisional frequency are of the same order of magnitude. The overall conclusion from these studies is that the single-fluid approach is suitable for dealing with slow processes in partially ionized plasmas, which happen on longer timescales than the collisional timescale. This approximation fails for timescales comparable to or shorter than the characteristic timescales for interaction between the charged and neutral particle populations.

Only recently, several numerical simulations of solar phenomena using a two-fluid approach have been reported. However, their realism is rather limited. Maneva et al. (2017) modeled magneto-acoustic wave propagation in the solar stratified chromosphere including effects of impact ionization and radiative recombination. These authors found numerous difficulties in constructing an equilibrium model atmosphere for the wave propagation, and in interpreting the results of their simulations in comparisons to more standard, single-fluid models. Magnetic reconnection has been studied using a two-fluid model by Smith & Sakai (2008), Leake et al. (2012, 2013), Alvarez Laguna et al. (2014), and Hillier et al. (2016). Important decoupling between the flows of neutrals and ions was observed in their simulations close to a reconnection site with considerable differences between inflows and outflows. Ni et al. (2018) also found that even in the absence of flow decoupling, dynamics of the ionization and recombination processes can have important consequences for the thermodynamics and observable properties of plasma at chromospheric reconnection sites. Martínez-Gomez et al. (2017) used a five-fluid model with three ionized and two neutral components, which takes into consideration Hall current and Ohm’s diffusion in addition to the friction due to collisions between different species. These latter authors apply their model to study the wave propagation in homogeneous plasmas composed of hydrogen and helium, similar to analytical models described above, but allowing time evolution of the ionization fraction. They confirmed the theoretical result that the inclusion of neutral components in the plasma modifies the oscillation periods of the low-frequency Alfvén waves and that collisions produce damping of the perturbations, the damping being more efficient when the collisional frequency is of the order of the oscillation frequency.

Considering all of the above it is evident that numerical modeling of solar chromospheric plasma including multi-fluid effects is a promising approach, and that a significant effort is required in this direction. Both previous theoretical and observational results suggest the need for such modeling. Here, we describe the results from the newly developed two-fluid code for purely hydrogen plasma which takes into account the effects of elastic and inelastic interactions between electrons, ions, and neutral atoms. In the current paper, which is the first in our study of chromospheric waves and shocks within a multi-fluid approximation, we verify the code against known analytical solutions for the wave propagation in a homogeneous plasma, and provide numerical convergence results of the implemented numerical scheme. In Sect. 5 we describe a test of a fast wave in a gravitationally stratified atmosphere, and test the numerical solution against the analytical solution. The full equations are provided in this first paper, describing the upgraded code, in order to enable reproducibility, and for future reference in the follow-on papers. In the second paper of the series (Popescu Braileanu et al. 2019), we apply the code to study the nonlinear behavior of chromospheric shocks under two-fluid conditions.

2. Numerical method

A multicomponent plasma evolution can be described by equations of continuity, momentum, and energy of its components which are obtained via momenta (0th, 1st, and 2nd moment) of the Boltzman equation. Typically, ions(i), electrons(e), and neutrals(n) are considered. When treated numerically with an explicit code, the equation of evolution of electrons imposes very small time steps and is therefore numerically restrictive. Because of the large ion-to-electron mass ratio, the electron propagation speed, sound speed, and Alfven velocity are also large relative to the corresponding bulk plasma parameters; this restricts the time step in an explicit scheme. However, for the same reason, the electron mass effects can generally be neglected from the electron equation of motion. This relaxes some of the numerical constraints and results in a generalized Ohm’s law used to compute the evolution of the electric field. The following assumptions are made in the set of two-fluid equations implemented in this work:

• Purely hydrogen plasma. This condition implies than only singly ionized ions are present and that charge neutrality is fulfilled on scales above the Debye radius scale, as usual. Charge exchange reactions are introduced via elastic collisions, by defining a charge-exchange collisional parameter α_cx, added to the elastic collisional parameter α. The expression of this term is given in Eq. (A.5). For such plasma, a constant mean molecular weight for neutrals is μ_n = 1 g mol⁻¹ and a constant mean molecular weight for charges is μ_c = 0.5 g mol⁻¹, where the sub-index “c” stands for “charges”.

• Identical temperatures of ions and electrons, that is, T_c = T_e = T_i. Neutrals and charges are allowed to have different temperatures, that is, generally T_n ≠ T_c.

• The center of mass velocity of charges is essentially ion velocity, that is, u_c = u_i. This condition implies that ${\mathit{u}}_{\mathrm{e}}={\mathit{u}}_{\mathrm{c}}-\frac{\mathit{J}}{e{n}_{\mathrm{e}}}$.

2.1. Equations

We solve continuity, momentum, and energy conservation equations written separately for charges (sub-index “c”) and neutrals (sub-index “n”). The interaction between neutrals and charges is introduced via the different kinds of collisional terms, as follows.

2.1.1. Continuity equations

$$\begin{array}{cc}\hfill \frac{\partial {\rho }_{\mathrm{n}}}{\partial t}+\mathrm{\nabla }·({\rho }_{\mathrm{n}}{\mathit{u}}_{\mathrm{n}})& =S_{\mathrm{n}},\hfill \\ \hfill \frac{\partial {\rho }_{\mathrm{c}}}{\partial t}+\mathrm{\nabla }·({\rho }_{\mathrm{c}}{\mathit{u}}_{\mathrm{c}})& =-S_{\mathrm{n}},\hfill \end{array}$$(1)

with:

$$\begin{array}{c}\hfill S_{\mathrm{n}}={\rho }_{\mathrm{c}}{\mathrm{\Gamma }}^{\mathrm{rec}}-{\rho }_{\mathrm{n}}{\mathrm{\Gamma }}^{\mathrm{ion}},\end{array}$$(2)

the latter being the ionization/recombination collisional source term. The variables Γ^ion and Γ^rec are ionization and recombination rates and are defined in Appendix A. The variables ρ_c and ρ_n are mass density of charges and neutrals, and u_c and u_n are their center of mass velocities.

2.1.2. Momentum equations

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

The elastic collisions between ions and neutrals and between electrons and neutrals can be combined into a single collision frequency between charges and neutrals. By defining the collisional parameter:

$$\begin{array}{c}\hfill \alpha =\frac{{\rho }_{\mathrm{e}}{\nu }_{\mathrm{en}}+{\rho }_{\mathrm{i}}{\nu }_{\mathrm{in}}}{{\rho }_{\mathrm{n}}{\rho }_{\mathrm{c}}},\end{array}$$(4)

the effective collision frequencies between charges and neutrals, and between neutrals and charges become, respectively:

$$\begin{array}{cc}\hfill {\nu }_{\mathrm{cn}}& =\alpha {\rho }_{\mathrm{n}},\hfill \\ \hfill {\nu }_{\mathrm{nc}}& =\alpha {\rho }_{\mathrm{c}}.\hfill \end{array}$$(5)

Here, ν_en and ν_in are the electron-neutral and ion-neutral collision frequencies, respectively, with expressions given in Eq. (A.3). If we take into account the charge-exchange reactions, α_cx, as defined in Eq. (A.5) should be added to the elastic collision parameter α. The momentum collisional term R_n can be written as

$$\begin{array}{c}\hfill {\mathit{R}}_{\mathrm{n}}={\rho }_{\mathrm{c}}{\mathit{u}}_{\mathrm{c}}{\mathrm{\Gamma }}^{\mathrm{rec}}-{\rho }_{\mathrm{n}}{\mathit{u}}_{\mathrm{n}}{\mathrm{\Gamma }}^{\mathrm{ion}}+\alpha {\rho }_{\mathrm{n}}{\rho }_{\mathrm{c}}({\mathit{u}}_{\mathrm{c}}-{\mathit{u}}_{\mathrm{n}}).\end{array}$$(6)

Here ${\widehat{\mathbf{p}}}_{\mathrm{c}}$ and ${\widehat{\mathbf{p}}}_{\mathrm{n}}$ are pressure tensors for charged and neutral species, B is magnetic field, J is current density, and R_n is momentum collisional exchange term. The latter takes into account momentum gains/losses due to elastic collisions, with collisional frequencies between ions, electrons, and neutrals, ν_in and ν_en, defined in Appendix A; it also takes into account momentum gains/losses due to ionization and recombination (the terms proportional to the ionization/recombination rates Γ).

We consider isotropic pressure, but we include viscosity in the pressure tensor:

$$\begin{array}{c}\hfill {\widehat{\mathbf{p}}}_{\alpha }=p_{\alpha }\mathbb{I}-{\widehat{\mathit{\tau }}}_{\alpha },\end{array}$$(7)

where the elements of the viscosity tensor are

$$\begin{array}{c}\hfill {{\tau }_{\alpha }}_{\mathit{ij}}={\xi }_{\alpha }(\frac{\partial {u_{\alpha }}_i}{\partial x_j}+\frac{\partial {u_{\alpha }}_j}{\partial x_i}).\end{array}$$(8)

The expressions used for the viscosity coefficient for species α, denoted by the symbol ξ_α, are given in Eq. (A.10).

The gravitational acceleration g is oriented in the negative z direction.

2.1.3. Energy conservation equations

$$\begin{array}{cc}& \frac{\partial }{\partial t}(e_{\mathrm{n}}+\frac{1}{2}{\rho }_{\mathrm{n}}{u}_{\mathrm{n}}^2)+\mathrm{\nabla }·({\mathit{u}}_{\mathrm{n}}(e_{\mathrm{n}}+\frac{1}{2}{\rho }_{\mathrm{n}}{u}_{\mathrm{n}}^2)+{\widehat{\mathbf{p}}}_{\mathrm{n}}·{\mathit{u}}_{\mathrm{n}}-K_{\mathrm{n}}\mathbf{\nabla }{T}_{\mathrm{n}})\hfill \\ \hfill & =-{\rho }_{\mathrm{n}}{\mathit{u}}_{\mathrm{n}}·\mathit{g}+M_{\mathrm{n}},\hfill \\ \hfill & \frac{\partial }{\partial t}(e_{\mathrm{c}}+\frac{1}{2}{\rho }_{\mathrm{c}}{u}_{\mathrm{c}}^2)+\mathrm{\nabla }·({\mathit{u}}_{\mathrm{c}}(e_{\mathrm{c}}+\frac{1}{2}{\rho }_{\mathrm{c}}{u}_{\mathrm{c}}^2)+{\widehat{\mathbf{p}}}_{\mathrm{c}}·{\mathit{u}}_{\mathrm{c}}-K_{\mathrm{c}}\mathbf{\nabla }{T}_{\mathrm{c}})\hfill \\ \hfill & =-{\rho }_{\mathrm{c}}{\mathit{u}}_{\mathrm{c}}·\mathit{g}+\mathit{J}·\mathit{E}-M_{\mathrm{n}},\hfill \end{array}$$(9)

with:

$$\begin{array}{cc}\hfill M_{\mathrm{n}}& =(\frac{1}{2}{\mathrm{\Gamma }}^{\mathrm{rec}}{\rho }_{\mathrm{c}}{u}_{\mathrm{c}}^2-\frac{1}{2}{\rho }_{\mathrm{n}}{u}_{\mathrm{n}}^2{\mathrm{\Gamma }}^{\mathrm{ion}})+\frac{1}{\gamma -1}\frac{k_{\mathrm{B}}}{m_{\mathrm{n}}}({\rho }_{\mathrm{c}}{T}_{\mathrm{c}}{\mathrm{\Gamma }}^{\mathrm{rec}}-{\rho }_{\mathrm{n}}{T}_{\mathrm{n}}{\mathrm{\Gamma }}^{\mathrm{ion}})\hfill \\ \hfill & +\frac{1}{2}({u_{\mathrm{c}}}^2-{u_{\mathrm{n}}}^2)\alpha {\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}+\frac{1}{\gamma -1}\frac{k_{\mathrm{B}}}{m_{\mathrm{n}}}(T_{\mathrm{c}}-T_{\mathrm{n}})\alpha {\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}.\hfill \end{array}$$(10)

In these equations, e_c and e_n are charges and neutral internal energy, respectively, and T_c and T_n are the corresponding temperatures, K_c and K_n are thermal conductivity coefficients, whose expressions are shown in Eq. (A.10). Among the four terms of the expression of M_n, the first two are related to ionization/recombination processes, and the last two terms are related to elastic collisions. Terms containing temperatures, namely the second and the fourth terms, represent the thermal exchange between charges and neutrals due to inelastic and elastic collisions, respectively. Because the energy equations evolve kinetic and internal energy, the M_n terms include kinetic energy gains/losses due to the change in the number of particles, and to the work done by the collisional terms. One can write equations of evolution of internal energy alone:

$$\begin{array}{cc}& \frac{\partial e_{\mathrm{n}}}{\partial t}+\mathrm{\nabla }·({\mathit{u}}_{\mathrm{n}}{e}_{\mathrm{n}}-K_{\mathrm{n}}\mathbf{\nabla }{T}_{\mathrm{n}})+p_{\mathrm{n}}\mathrm{\nabla }·{\mathit{u}}_{\mathrm{n}}-{\widehat{\tau }}_{\mathrm{n}}:\mathbf{\nabla }{\mathit{u}}_{\mathrm{n}}=M_{\mathrm{n}}^{\prime },\hfill \\ \hfill & \frac{\partial e_{\mathrm{c}}}{\partial t}+\mathrm{\nabla }·({\mathit{u}}_{\mathrm{c}}{e}_{\mathrm{c}}-K_{\mathrm{c}}\mathbf{\nabla }{T}_{\mathrm{c}})+p_{\mathrm{c}}\mathrm{\nabla }·{\mathit{u}}_{\mathrm{c}}-{\widehat{\tau }}_{\mathrm{c}}:\mathbf{\nabla }{\mathit{u}}_{\mathrm{c}}=\mathit{J}·{\mathit{E}}_{\mathrm{diff}}+M_{\mathrm{c}}^{\prime },\hfill \end{array}$$(11)

where

$$\begin{array}{c}\hfill {\mathit{E}}_{\mathrm{diff}}=\mathit{E}+{\mathit{u}}_{\mathrm{c}}\times \mathit{B}\end{array}$$

is the nonideal part of the electric field E, and

$$\begin{array}{c}\hfill {\widehat{\tau }}_{\alpha }:\mathbf{\nabla }\mathit{u}={\xi }_{\alpha }{\displaystyle \sum _{i<j}}{(\frac{\partial {u_{\alpha }}_i}{\partial x_j}+\frac{\partial {u_{\alpha }}_j}{\partial x_i})}^2,\end{array}$$

which is the heating due to the viscosity, is a non-negative quantity. Subsequently, the collisional terms become, for neutrals and charges:

$$\begin{array}{cc}\hfill M_{\mathrm{n}}^{\prime }=& M_{\mathrm{n}}-{\mathit{u}}_{\mathrm{n}}·{\mathit{R}}_{\mathrm{n}}+\frac{1}{2}{u}_{\mathrm{n}}^2{S}_{\mathrm{n}}=\frac{1}{2}{\mathrm{\Gamma }}^{\mathrm{rec}}{\rho }_{\mathrm{c}}{({\mathit{u}}_{\mathrm{c}}-{\mathit{u}}_{\mathrm{n}})}^2\hfill \\ \hfill & +\frac{1}{\gamma -1}\frac{k_{\mathrm{B}}}{m_{\mathrm{n}}}({\rho }_{\mathrm{c}}{T}_{\mathrm{c}}{\mathrm{\Gamma }}^{\mathrm{rec}}-{\rho }_{\mathrm{n}}{T}_{\mathrm{n}}{\mathrm{\Gamma }}^{\mathrm{ion}})\hfill \\ \hfill & +\frac{1}{2}{({\mathit{u}}_{\mathrm{c}}-{\mathit{u}}_{\mathrm{n}})}^2\alpha {\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}\hfill \\ \hfill & +\frac{1}{\gamma -1}\frac{k_{\mathrm{B}}}{m_{\mathrm{n}}}(T_{\mathrm{c}}-T_{\mathrm{n}})\alpha {\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}},\hfill \\ \hfill M_{\mathrm{c}}^{\prime }=& -M_{\mathrm{n}}+{\mathit{u}}_{\mathrm{c}}·{\mathit{R}}_{\mathrm{n}}-\frac{1}{2}{u}_{\mathrm{c}}^2{S}_{\mathrm{n}}=\frac{1}{2}{\mathrm{\Gamma }}^{\mathrm{ion}}{\rho }_{\mathrm{n}}{({\mathit{u}}_{\mathrm{c}}-{\mathit{u}}_{\mathrm{n}})}^2\hfill \\ \hfill & -\frac{1}{\gamma -1}\frac{k_{\mathrm{B}}}{m_{\mathrm{n}}}({\rho }_{\mathrm{c}}{T}_{\mathrm{c}}{\mathrm{\Gamma }}^{\mathrm{rec}}-{\rho }_{\mathrm{n}}{T}_{\mathrm{n}}{\mathrm{\Gamma }}^{\mathrm{ion}})\hfill \\ \hfill & +\frac{1}{2}{({\mathit{u}}_{\mathrm{c}}-{\mathit{u}}_{\mathrm{n}})}^2\alpha {\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}\hfill \\ \hfill & -\frac{1}{\gamma -1}\frac{k_{\mathrm{B}}}{m_{\mathrm{n}}}(T_{\mathrm{c}}-T_{\mathrm{n}})\alpha {\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}.\hfill \end{array}$$(12)

We observe that the second and the fourth terms remain unchanged. However, the first and the third terms are different; they become non-negative quantities and thus contribute to dissipative heating. The neutrals and charges are heated because of recombination and ionization processes, respectively (the first term), and both species are heated at the same rate by elastic collisions (the third term, equal for neutrals and charges, also called “frictional heating”). Most expressions of the terms related to collisions (including coefficients of viscosities, thermal conductions, and magnetic resistivity) are derived by Braginskii (1965), but they are also detailed in the descriptions of models by Meier & Shumlak (2012), Leake et al. (2012, 2013).

The relation between internal energy and pressure, and between the pressure, the number density, and the temperature are defined by ideal gas laws:

$$\begin{array}{cc}\hfill e_{\mathrm{c},\mathrm{n}}& =p_{\mathrm{c},\mathrm{n}}/(\gamma -1),\hfill \\ \hfill p_{\mathrm{c},\mathrm{n}}& =n_{\mathrm{c},\mathrm{n}}{k}_{\mathrm{B}}{T}_{\mathrm{c},\mathrm{n}}.\hfill \end{array}$$(13)

2.1.4. Ohm’s law

In order to obtain Ohm’s law we operate on the momentum equation for ions and electrons. For the case of purely hydrogen plasma, and by neglecting terms proportional to m_e/m_i, such as electron inertia terms, and time variation of currents (see Khomenko et al. 2014; Ballester et al. 2018), we obtain to the following expression.

$$\begin{array}{c}\hfill \mathit{E}+{\mathit{u}}_{\mathrm{c}}\times \mathit{B}={\eta }_{\mathrm{H}}\mathit{J}\times \mathit{B}-{\eta }_{\mathrm{H}}\mathbf{\nabla }{p}_{\mathrm{e}}+\eta \mathit{J}-{\eta }_{\mathrm{D}}({\mathit{u}}_{\mathrm{c}}-{\mathit{u}}_{\mathrm{n}}),\end{array}$$(14)

where

$$\begin{array}{cc}& {\eta }_{\mathrm{H}}=\frac{1}{e{n}_{\mathrm{e}}},\hfill \\ \hfill & \eta =\frac{{\rho }_{\mathrm{e}}({\nu }_{\mathrm{en}}+{\nu }_{\mathrm{ei}})}{{(e{n}_{\mathrm{e}})}^2},\hfill \\ \hfill & {\eta }_{\mathrm{D}}=\frac{{\rho }_{\mathrm{e}}({\nu }_{\mathrm{en}}-{\nu }_{\mathrm{in}})}{e{n}_{\mathrm{e}}}.\hfill \end{array}$$(15)

All four terms on the right-hand side of Eq. (14) are the same as in the single fluid approach, and the first three of them represent the Hall term, the Biermann battery term, and the Ohmic term. The fourth term, proportional to the difference of the velocities of the two species, is usually neglected in both approximations. The ambipolar term, which appears in the single-fluid approximation as a consequence of the charged and neutral fluids having different center-of-mass velocities does not appear in the two-fluid approach. The derivation of these terms can be found in Khomenko et al. (2014).

The induction equation for the evolution of the magnetic field is obtained using Maxwell equations, in a usual way,

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}=-\mathbf{\nabla }\times \mathit{E}.\end{array}$$(16)

2.2. Numerical implementation

2.2.1. MANCHA3D code

In this work we have extended the MANCHA3D code, a single-fluid code, where the effects of partial ionization have been implemented through a generalized Ohm’s law including Ohmic, ambipolar, and Hall and Biermann battery terms. The new code solves the two-fluid equations in approximation of purely hydrogen plasma, as described above. Ions and electrons are evolved together through the same continuity, momentum, and energy equations, and neutrals are evolved separately. The neutral and charges equations are coupled by collisions. These new developments are described below.

The original MANCHA3D code is described elsewhere (Khomenko & Collados 2006, 2012; Felipe et al. 2010; González-Morales et al. 2018). The code is fully 3D, written in Fortran 90, and parallelized with MPI; the output files are in HDF5 format. The units used in the equations are SI units. The variables are split into background and perturbation variables, keeping all nonlinear terms. The viscous forces and the Ohmic dissipation tend to eliminate high-frequency oscillations. However, collisional viscosity (ξ_c, n; Eq. (8)) and Ohmic (ν; Eq. (14)) coefficients for the parameters of the solar atmosphere are too small to effectively damp nonresolved variations of different parameters at the grid cell size. Therefore, for numerical stability, these coefficients have corresponding numerical analogs referred to as hyperdiffusive coefficients. The hyperdiffusive coefficients correspond to physical terms: magnetic Ohmic diffusion, conductivity, viscosity, and so on. The amplitude of the hyperdiffusive coefficients is a complex function of space coordinates and of gradients of the corresponding variables (for details see Vögler et al. 2005; Felipe et al. 2010). Since high numerical diffusion is not always desirable, an additional filtering of small wavelengths following Parchevsky & Kosovichev (2007) is implemented in MANCHA3D. In order to avoid reflections at the boundaries, MANCHA3D uses the Perfect Matching Layer technique (Felipe et al. 2010). The spatial discretization scheme is fourth-order accurate. The original temporal discretization scheme was a fourth-order explicit Runge–Kutta, modified to an arbitrary order in the latest version of the code due to implementation of STS and HDS schemes (González-Morales et al. 2018). The two-fluid implementation required further modification of the explicit Runge–Kutta integration. The new implementation is a semi-implicit scheme with up to second-order accuracy, as is shown in the numerical tests of the code described below.

2.2.2. Numerical scheme

Since the solar atmosphere is strongly gravitationally stratified, the collisional frequencies change orders of magnitude from the photosphere to the chromosphere. In the deep layers, collisional terms are typically very large and dominate over the hydrodynamical timescales. In the chromosphere this situation may change, since collisional frequencies become significantly smaller. Therefore, the numerical scheme of our two-fluid numerical code has to be able to handle large variations in collision frequencies. In the case of dominant collisional terms the equations become stiff and therefore an explicit integration scheme is unstable. This is because the collisional terms that appear in the equations of continuity, momentum, and energy are proportional to the collision frequency, and the time step should be at most the inverse of the collision frequency in order to ensure stability in an explicit scheme. This imposes very small time steps when the collisions are dominant. A suitable alternative is to use semi-implicit schemes similar to that described by Toth et al. (2012), such as the one implemented here. Our system of equations, Eqs. (1), (3), (9), and (14), can be schematically written as follows.

$$\begin{array}{c}\hfill \frac{\partial \mathit{U}}{\partial t}=\mathit{R}(\mathit{U}),\end{array}$$(17)

where U is the vector of variables,

$$\begin{array}{c}\hfill \mathit{U}=\{{\rho }_{\mathrm{n}},{\rho }_{\mathrm{c}},{\rho }_{\mathrm{n}}{\mathit{u}}_{\mathrm{n}},{\rho }_{\mathrm{c}}{\mathit{u}}_{\mathrm{c}},e_{\mathrm{n}}+\frac{1}{2}{\rho }_{\mathrm{n}}{u}_{\mathrm{n}}^2,e_{\mathrm{c}}+\frac{1}{2}{\rho }_{\mathrm{c}}{u}_{\mathrm{c}}^2,\mathit{B}\},\end{array}$$(18)

and R(U) is the sum of fluxes and collisional terms.

In order to evolve this system of equations in time: we split R into an explicit operator E and an implicit operator P,

$$\begin{array}{c}\hfill \mathit{R}(\mathit{U})=\mathit{E}(\mathit{U})+\mathit{P}(\mathit{U}).\end{array}$$(19)

The explicit operator E(U) contains fluxes, while the implicit operator contains the collisional terms, Eqs. (2), (6), and (10). With U^* being the solution after applying the explicit update,

$$\begin{array}{c}\hfill {\mathit{U}}^{\ast }={\mathit{U}}^n+\mathrm{\Delta }t\mathit{E},\end{array}$$(20)

we apply the implicit update using the following discretization.

$$\begin{array}{c}\hfill {\mathit{U}}^{n+1}={\mathit{U}}^{\ast }+(1-\beta )\mathrm{\Delta }t\mathit{P}({\mathit{U}}^{\ast })+\beta \mathrm{\Delta }t\mathit{P}({\mathit{U}}^{n+1}).\end{array}$$(21)

The value of β can vary between 0 and 1, and its influence on the final solution will be discussed below.

The code solves three systems of two coupled equations during the implicit update, with U being in the three cases a vector with two components, as defined above by Eq. (18), $\mathit{U}=\{\{{\rho }_{\mathrm{n}},{\rho }_{\mathrm{c}}\},\{{\rho }_{\mathrm{n}}{\mathit{u}}_{\mathrm{n}},{\rho }_{\mathrm{c}}{\mathit{u}}_{\mathrm{c}}\},\{e_{\mathrm{n}}+\frac{1}{2}{\rho }_{\mathrm{n}}{u}_{\mathrm{n}}^2,e_{\mathrm{c}}+\frac{1}{2}{\rho }_{\mathrm{c}}{u}_{\mathrm{c}}^2\}\}$, meaning that all of the collisional terms, namely S_n, R_n, and M_n defined in Eqs. (2), (6), and (10), can be treated in the same way. We note however that depending on the magnitude of the collisional terms S_n in the continuity equations, it might not be worth treating them implicitly. In such cases we add these terms into the explicit operator E(U). All the terms on the right-hand side of the induction equation, Eq. (14), are treated explicitly.

We linearize P(U^n + 1) around the value of U^* after an explicit update:

$$\begin{array}{c}\hfill \mathit{P}({\mathit{U}}^{n+1})=\mathit{P}({\mathit{U}}^{\ast })+\widehat{J}·({\mathit{U}}^{n+1}-{\mathit{U}}^{\ast }),\end{array}$$

where the Jacobian, Ĵ, is defined as follows.

$$\begin{array}{c}\hfill \widehat{J}=\frac{\partial \mathit{P}}{\partial \mathit{U}}·\end{array}$$(22)

Due to particular symmetry properties of the equations, J₂₁ = −J₁₁, J₂₂ = −J₁₂. This allows the solution for the implicit update to be calculated directly. The update in one sub-step of a two-step Runge–Kutta scheme can then be summarized as

$$\begin{array}{c}\hfill {\mathit{U}}^{n+\frac{1}{k}}={\mathit{U}}^{\ast }+\frac{\mathrm{\Delta }{t}_{\frac{1}{k}}((1-{\beta }_{\frac{1}{k}}){\mathit{P}}^n+{\beta }_{\frac{1}{k}}{\mathit{P}}^{\ast })}{1+{\beta }_{\frac{1}{k}}\mathrm{\Delta }{t}_{\frac{1}{k}}({\widehat{J}}_{12}-{\widehat{J}}_{11})},\end{array}$$(23)

where k = 2, 1; ${\mathrm{\Delta }}_{\frac{1}{k}}=\frac{1}{k}\mathrm{\Delta }t$ is the time step used in the substep, and Δt is the complete time step. The scheme parameter ${\beta }_{\frac{1}{k}}$ can have a different value in each substep (subscript $\frac{1}{k}$). The magnitude of Δt is limited by the MHD CFL condition. The explicit operator, E, is recalculated in each Runge–Kutta sub-step using the variables updated after the previous sub-step. The Jacobian Ĵ is also calculated after the explicit update, so that Ĵ = Ĵ^∗.

The elements of the Jacobian can be written in the analytical form.

For the momentum equations, this ensures that

$$\begin{array}{c}\hfill {\mathit{P}}^{\ast }=\widehat{J}·{\mathit{U}}^{\ast }.\end{array}$$(24)

For the energy equations, the same expression becomes

$$\begin{array}{c}\hfill {\mathit{P}}^{\ast }=\widehat{J}·{\mathit{U}}_{\mathrm{T}}^{\ast ,1}+\widehat{K}·{\mathit{U}}_{\mathrm{m}}^{n+\frac{1}{k}},\end{array}$$(25)

where U_m and U_T are the kinetic and internal energy parts of the variable U, respectively. The Ĵ and K̂ for the energy equation are defined as follows.

$$\begin{array}{c}\hfill \widehat{J}=\frac{\partial \mathit{P}}{\partial {\mathit{U}}_{\mathrm{T}}};\widehat{K}=\frac{\partial \mathit{P}}{\partial {\mathit{U}}_{\mathrm{m}}}·\end{array}$$(26)

The variable ${\mathit{U}}_{\mathrm{T}}^{\ast ,1}$ is the difference between U^* after an explicit update, and one Runge–Kutta sub-step of the implicit update, i.e.,

$$\begin{array}{c}\hfill {\mathit{U}}_{\mathrm{T}}^{\ast ,1}={\mathit{U}}^{\ast }-{\mathit{U}}_{\mathrm{m}}^{n+\frac{1}{k}}.\end{array}$$(27)

Overall, this means that the energy implicit update is done with variables calculated after the momentum implicit update. The analytical expressions for Ĵ and K̂ are given in Appendix B. The numerical scheme is stable when we use the time step imposed by the CFL. For ${\beta }_{\frac{1}{2}}=1$ and ${\beta }_1=\frac{1}{2}$ the scheme is formally second-order accurate.

3. Analytical model of waves in a uniform plasma

Our newly implemented scheme has been checked using simulations of wave propagation in a uniform plasma. We used two simple analytical cases. In the first case we consider purely acoustic waves propagating in an atmosphere where the background equilibrium temperature of charges and neutrals is different. In this case the difference in behavior between charges and neutrals is created by the different sound speeds of the background atmosphere. In the second case we consider the propagation of Alfvén waves. Here the difference in the behavior of charges and neutrals is produced additionally due to the presence of the magnetic field. In this simple setup the problems have analytical solutions that can be obtained by linearizing the two fluid equations and looking for solutions in the form of monochromatic waves. The sections below describe the analytical solution and results for both sets of simulations, together with the tests of numerical convergence of the scheme.

3.1. Acoustic waves

Acoustic waves are longitudinal waves, therefore their velocity vector and direction of propagation are parallel, that is, u_n ∥ k and u_c ∥ k. We choose the direction of propagation along z direction. We neglect the effects of ionization, recombination, thermal exchange, thermal conduction, and radiation. Therefore, we only consider the momentum exchange due to elastic collisions. In this simple situation the linearized continuity, momentum, and adiabatic energy equations become

$$\begin{array}{cc}& \frac{\partial {\rho }_{\mathrm{n}1}}{\partial t}+{\rho }_{\mathrm{n}0}\frac{\partial u_{\mathrm{nz}}}{\partial z}=0,\hfill \\ \hfill & \frac{\partial {\rho }_{\mathrm{c}1}}{\partial t}+{\rho }_{\mathrm{c}0}\frac{\partial u_{\mathrm{cz}}}{\partial z}=0,\hfill \\ \hfill & {\rho }_{\mathrm{n}0}\frac{\partial u_{\mathrm{nz}}}{\partial t}+\frac{\partial p_{\mathrm{n}1}}{\partial z}=\alpha {\rho }_{\mathrm{n}0}{\rho }_{\mathrm{c}0}(u_{\mathrm{cz}}-u_{\mathrm{nz}}),\hfill \\ \hfill & {\rho }_{\mathrm{c}0}\frac{\partial u_{\mathrm{cz}}}{\partial t}+\frac{\partial p_{\mathrm{c}1}}{\partial z}=-\alpha {\rho }_{\mathrm{n}0}{\rho }_{\mathrm{c}0}(u_{\mathrm{cz}}-u_{\mathrm{nz}}),\hfill \\ \hfill & \frac{\partial p_{\mathrm{n}1}}{\partial t}-{c_{\mathrm{n}0}}^2\frac{\partial {\rho }_{\mathrm{n}1}}{\partial t}=0,\hfill \\ \hfill & \frac{\partial p_{\mathrm{c}1}}{\partial t}-{c_{\mathrm{c}0}}^2\frac{\partial {\rho }_{\mathrm{c}1}}{\partial t}=0,\hfill \end{array}$$(28)

where u_nz and u_cz are the velocity projections along z direction, and c_n0 and c_c0 are sound speeds of neutrals and charges. Variables with subscript “0” refer to the background unperturbed atmosphere, and variables with subscript “1” refer to the perturbation. The collisional parameter α defined by Eq. (4) uses the values of the unperturbed mass densities of electrons, ions, neutrals, and charges. Since for the simple model assumed here the value of α only depends on the homogeneous background variables in the linear approximation, it is kept constant in time and space.

We search for solutions of form,

$$\begin{array}{c}\hfill \{u_{\mathrm{nz}},u_{\mathrm{cz}},{\rho }_{\mathrm{n}1},{\rho }_{\mathrm{c}1},p_{\mathrm{n}1},p_{\mathrm{c}1}\}=\{U_{\mathrm{n}},U_{\mathrm{c}},R_{\mathrm{n}},R_{\mathrm{c}},P_{\mathrm{n}},P_{\mathrm{c}}\}{\mathrm{e}}^{i(\omega t-kz)},\end{array}$$(29)

where {U_n, U_c, R_n, R_c, P_n, P_c} are, in general, complex amplitudes for all the perturbed variables. These amplitudes are related through the so-called polarization relation, so one needs to set one of the amplitudes to obtain the rest of them through the following relations:

$$\begin{array}{cc}\hfill U_{\mathrm{c}}& =\frac{\omega R_{\mathrm{c}}}{k{\rho }_{\mathrm{c}0}};U_{\mathrm{n}}=\frac{\omega R_{\mathrm{n}}}{k{\rho }_{\mathrm{n}0}},\hfill \\ \hfill \frac{R_{\mathrm{c}}}{R_{\mathrm{n}}}& =\frac{i{\omega }^2-i{k}^2{c}_{\mathrm{n}0}^2+\alpha {\rho }_{\mathrm{c}0}\omega }{\alpha {\rho }_{\mathrm{n}0}\omega },\hfill \\ \hfill P_{\mathrm{c}}& =c_{\mathrm{c}0}^2{R}_{\mathrm{c}};P_{\mathrm{n}}=c_{\mathrm{n}0}^2{R}_{\mathrm{n}}.\hfill \end{array}$$(30)

The resulting dispersion relation, which relates the frequency and the wavenumber, has the following form:

$$\begin{array}{cc}\hfill i{\omega }^4+\alpha {\omega }^3({\rho }_{\mathrm{c}0}+{\rho }_{\mathrm{n}0})& -i{\omega }^2{k}^2(c_{\mathrm{n}0}^2+c_{\mathrm{c}0}^2)-\omega k^2\alpha ({\rho }_{\mathrm{c}0}{c}_{\mathrm{c}0}^2\hfill \\ \hfill & +{\rho }_{\mathrm{n}0}{c}_{\mathrm{n}0}^2)+i{k}^4{c}_{\mathrm{n}0}^2{c}_{\mathrm{c}0}^2=0.\hfill \end{array}$$(31)

Once the background values are specified, the value of the collisional parameter α is determined. It is then convenient to operate in terms of the following parameters: ω/k, and α/k. The latter ratio is an effective measure of the collisional strength. The analysis below will be done in terms of the nondimensional variables:

$$\begin{array}{c}\hfill E=\frac{\omega }{k{c}_{\mathrm{tot}}};F=\frac{\alpha {\rho }_{\mathrm{tot}}}{k{c}_{\mathrm{tot}}}·\end{array}$$(32)

These variables are similar to those used in Soler et al. (2013a). Later on in the paper we represent the solution of the dispersion relation in terms of these nondimensional variables, keeping in mind that F is inversely proportional to the wavenumber k with the mutiplicative factors fixed (by “varying the collisional strength F” we understand varying the wavenumber k, since the value of α is fixed).

Using nondimensional variables, the polarization relations and the dispersion relation become:

$$\begin{array}{c}\hfill \frac{U_{\mathrm{c}}}{R_{\mathrm{c}}}=E\frac{c_{\mathrm{tot}}}{{\rho }_{\mathrm{c}0}};\frac{R_{\mathrm{n}}}{R_{\mathrm{c}}}=\frac{EF{\rho }_{\mathrm{n}0}/{\rho }_{\mathrm{tot}}}{i{E}^2-i{c_{\mathrm{n}}}_0^2/c_{\mathrm{tot}}^2+EF{\rho }_{\mathrm{c}0}/{\rho }_{\mathrm{tot}}},\end{array}$$(33)

$$\begin{array}{c}\hfill i{E}^4+F{E}^3-i{E}^2{a}_2-EF+i{a}_0=0,\end{array}$$(34)

where:

$$\begin{array}{c}\hfill a_2=\frac{c_{\mathrm{n}0}^2+c_{\mathrm{c}0}^2}{c_{\mathrm{tot}}^2};a_0=\frac{c_{\mathrm{n}0}^2{c}_{\mathrm{c}0}^2}{c_{\mathrm{tot}}^4}·\end{array}$$(35)

Here ρ_tot = ρ_n0 + ρ_c0 and $c_{\mathrm{tot}}=\sqrt{\gamma (p_{\mathrm{n}0}+p_{\mathrm{c}0})/({\rho }_{\mathrm{n}0}+{\rho }_{\mathrm{c}0})}$ are the total density and the sound speed of the whole fluid.

3.2. Alfvén waves

Alfvén waves are transverse incompressible waves, polarized perpendicularly to the direction of the background magnetic field. The analytical solution for the Alfvén wave propagation in partially ionized plasma used in our work is similar to those developed previously by Soler et al. (2013a). We consider here a simple case of Alfvén waves with propagation only along the background magnetic field. These are frequently called pure Alfvén waves. In our case, k ∥ B₀, u_n ∥ B₁, u_c ∥ B₁, B₁ ⊥ B₀. We choose the equilibrium magnetic field along z direction and perturbation of magnetic field and velocities along x direction.

$$\begin{array}{cc}\hfill {\mathit{B}}_0& =(0,0,B_0),\mathit{k}=(0,0,k),{\mathit{B}}_1=(B_{1x},0,0),\hfill \\ \hfill {\mathit{u}}_c& =(u_{\mathrm{cx}},0,0),{\mathit{u}}_n=(u_{\mathrm{nx}},0,0).\hfill \end{array}$$

Otherwise the assumptions of the equations are similar to the acoustic wave case. The linearized momentum and induction equations take the following form.

$$\begin{array}{cc}& {\rho }_{\mathrm{n}0}\frac{\partial u_{\mathrm{nx}}}{\partial t}=\alpha {\rho }_{\mathrm{n}0}{\rho }_{\mathrm{n}0}(u_{\mathrm{cx}}-u_{\mathrm{nx}}),\hfill \\ \hfill & {\rho }_{\mathrm{c}0}\frac{\partial u_{\mathrm{cx}}}{\partial t}=\frac{B_0}{{\mu }_0}\frac{\partial B_{1x}}{\partial z}-\alpha {\rho }_{\mathrm{c}0}{\rho }_{\mathrm{n}0}(u_{\mathrm{cx}}-u_{\mathrm{nx}}),\hfill \\ \hfill & \frac{\partial B_{1x}}{\partial t}=B_0\frac{\partial u_{\mathrm{cx}}}{\partial z}·\hfill \end{array}$$(36)

The equations of continuity and energy conservation are not needed in this case because the Alfvén waves are purely incompressible. Similar to acoustic waves we use the solutions of the form

$$\begin{array}{c}\hfill \{u_{\mathrm{nx}},u_{\mathrm{cx}},B_{1x}\}=\{U_{\mathrm{n}},U_{\mathrm{c}},B\}{\mathrm{e}}^{i(\omega t-kz)},\end{array}$$(37)

where the amplitudes can be complex in general. Plugging this ansatz into Eq. (36), we obtain the following polarization relations.

$$\begin{array}{cc}& B=-\frac{k{U}_{\mathrm{c}}{B}_0}{\omega },\hfill \\ \hfill & \frac{U_{\mathrm{c}}}{U_{\mathrm{n}}}=\frac{i{\rho }_{\mathrm{n}0}\omega +\alpha {\rho }_{\mathrm{n}0}{\rho }_{\mathrm{c}0}}{\alpha {\rho }_{\mathrm{n}0}{\rho }_{\mathrm{c}0}}·\hfill \end{array}$$(38)

and the dispersion relation in the following form.

$$\begin{array}{c}\hfill -{\mu }_0{\rho }_{\mathrm{c}0}{\omega }^3+i{\mu }_0\alpha {\rho }_{\mathrm{c}0}({\rho }_{\mathrm{c}0}+{\rho }_{\mathrm{n}0}){\omega }^2+{B_0}^2{k}^2\omega -i\alpha {\rho }_{\mathrm{c}0}{B_0}^2{k}^2=0.\end{array}$$(39)

For the same reason as for the acoustic waves, we vary the collision strength by varying the wavenumber k parameter to study its influence on the propagation properties of the waves. Using the nondimensional variables:

$$\begin{array}{c}\hfill C=\frac{\omega }{k{v_{\mathrm{A}}}_0};D=\frac{\alpha {\rho }_{\mathrm{tot}}}{k{v_{\mathrm{A}}}_0},\end{array}$$(40)

the dispersion relation, and the polarization relations are

$$\begin{array}{c}\hfill -C^3+iD{C}^2+C-iD\frac{{\rho }_{\mathrm{c}0}}{{\rho }_{\mathrm{tot}}}=0,\end{array}$$(41)

$$\begin{array}{c}\hfill \frac{B}{U_{\mathrm{c}}}=-\frac{B_0}{C{v_{\mathrm{A}}}_0},\frac{U_{\mathrm{c}}}{U_{\mathrm{n}}}=1+i\frac{C}{D}\frac{{\rho }_{\mathrm{tot}}}{{\rho }_{\mathrm{c}0}}·\end{array}$$(42)

3.3. Background atmosphere

The values of the uniform background atmosphere used in this study are given in Table 1. For the Alfvén wave tests, a z-directed magnetic field B₀ = 5 × 10⁻³T is also added. The corresponding Alfvén speed calculated using the density of charges is ${{\mathit{v}}_{\mathrm{A}}}_0=B_0/\sqrt{{\mu }_0{\rho }_{\mathrm{c}0}}=6.32\times {10}^4$ m s⁻¹. The collisional parameter α defined in Eq. (4) calculated for the background has the value 7.324 × 10¹² m³ kg⁻¹ s⁻¹. In this test the temperature of electrons is assumed to be zero, they do not contribute to the pressure of charges, and the effective neutral-charge and charge-neutral collision frequencies defined in Eq. (5) become the neutral-ion and ion-neutral collision frequencies, with values ν_ni ≈ 3.65 × 10⁴ s⁻¹, and ν_in ≈ 7.3 × 10⁴ s⁻¹, respectively.

3.4. Parameters of the perturbation

We have to supply several free parameters to calculate the solution. We have chosen the real wave number k, and obtain the complex ω from the dispersion relations, Eqs. (31) and (39). This means that we consider wave damping in time. Additionally, since the perturbations in all the magnitudes are related, we need to supply the amplitude of one of them in order to calculate the others. Generally, since the polarization relations are complex, there will be the phase shift between oscillations of different quantities.

We have selected the wavenumber to be k = 2πn/L_z m⁻¹, with n being an integer number. We vary the wavenumber k by varying the domain length L_z in order to have n = 11 (a random choice) for all the simulations.

For the acoustic waves we choose the amplitude of the perturbation in density of charges, R_c = 10⁻³ρ_c0, and for the Alfvén waves we choose the amplitude of the perturbation in velocity of charges, U_c = 10⁻⁵v_A0. We calculate the other amplitudes from Eq. (30) for the acoustic waves and Eq. (36) for the Alfvén waves. In all cases, for the numerical solution we used periodic boundary conditions. For the initial condition of the perturbation we used the analytical solution at t = 0 s in the whole domain.

We run all simulations for a total time t_F so that we have L_z/t_F equal to 7447.273 m s⁻¹ for all the simulations of acoustic waves and equal to 32125.49 m s⁻¹ for the Alfvén wave simulations. These values are of the order of the characteristic speeds: the sound speed and the Alfvén speed.

3.5. Solution of the dispersion relation for acoustic waves

The acoustic wave dispersion relation is a fourth-order equation in both ω and k. If the neutrals and charges have the same background temperature, the dispersion relation has four solutions, with simple expressions as shown in Eq. (44). In our case T_c0/T_n0 = 4, and the values of the coefficients in Eq. (34) are a₂ = 5/2 and a₀ = 1, and the expressions for the solutions, obtained using Mathematica software, in this case are as follows.

$$\begin{array}{cc}\hfill E_{4,1}& =\frac{1}{4}(iF-\sqrt{2-F^2}\pm \sqrt{18-2F^2-2iF\sqrt{2-F^2}}),\hfill \\ \hfill E_{3,2}& =\frac{1}{4}(iF+\sqrt{2-F^2}\pm \sqrt{18-2F^2+2iF\sqrt{2-F^2}}).\hfill \end{array}$$(43)

We compute the complex frequency, ω, for the values of k between 6180 and 6.18 × 10⁻⁶ m⁻¹, for each of the four expressions of the mathematical solutions shown in Eq. (43), and marked as solutions 1, 2, 3, and 4 in Fig. 1. The real and imaginary parts of the four solutions are shown in the normalized units, E = f(F), defined in Eq. (32), that is, ω_{I, R}/k/c_tot as a function of αρ_tot/k/c_tot for better visualization. The solutions with positive ω_R travel in the positive direction of the z axis, and those with negative ω_R travel in the negative direction.

Fig. 1. Solutions of the dispersion relation, Eq. (31), for the acoustic wave case. Left and right panels present the real and imaginary parts of the dimensionless wave frequency ω, respectively, as functions of the dimensionless wavenumber k. The axes are scaled in nondimensional units E and F, defined in Eq. (32) for better visualization. Four solutions are marked with different colors and line styles: solid red for solution 1, orange dashed for solution 2, blue dot-dashed for solution 3, and green dotted for solution 4. We note that solutions 1 and 4 have the same imaginary parts and are superposed in the right-hand panel. Solutions 2 and 3 have the same real part equal to zero for large α/k and are also superimposed on top of one another. The values of k for which we compare the analytical and numerical solution using the values of ω, corresponding to solution 4 (green dots), are marked in the panels: orange circle (k = 618 m−1), pink triangle (k = 6.18 m−1), cyan diamond (k = 2.06 m−1), and black “X” (k = 6.18 × 10−6 m−1).

Waves corresponding to solutions 1 and 4 (red and green colors) propagate with a speed (v_ph = ω_R/k) equal to the sound speed of the charges ($c_{\mathrm{c}0}/c_{\mathrm{tot}}=\sqrt{2}\approx 1.414$) for small ratios α/k, in the negative and positive direction of the z axis, respectively. For the large values of α/k, their propagation speed becomes that of the whole fluid, c_tot. The imaginary part of the frequency ω_I is positive meaning wave damping. The imaginary part would be zero if the charges and neutrals had the same background temperature, since in that case the solutions of the dispersion relation, Eq. (34), would simplify to

$$\begin{array}{c}\hfill E_{4,1}=\pm 1,\end{array}$$(44)

with zero imaginary part for solutions 1 and 4. As the only difference between neutrals and charges in this particular model is in their sound speed, it is not surprising that by making the sound speed the same the damping disappears, since both fluids would oscillate with exactly the same velocity and therefore the frictional damping term would become strictly zero.

The value of damping is the same for solutions 1 and 4. For either very low or very high values of the ratio α/k, the damping relative to the wavenumber (ω_I/k) is small and approaching zero. The value of the ratio ω_I/k is maximum at a point located on the x axis at $F_{\mathrm{M}}=\alpha {\rho }_{\mathrm{tot}}/k_{\mathrm{M}}/c_{\mathrm{tot}}=\sqrt{2}$, and corresponds to k_M ≈ 7.5 m⁻¹. E(F = F_M) = $1/4(\sqrt{14}+i\sqrt{2})$; this gives the value of ${\omega }_{\mathrm{R}}(F=F_{\mathrm{M}})/(\alpha {\rho }_{\mathrm{n}0})=3\sqrt{7}/8\approx 1$, which means that the real part of the wave frequency is approximately equal to the ion-neutral collision frequency at the point where the damping relative to wavenumber, the ratio ω_I/k, is maximum (F_M).

Waves corresponding to solutions 2 and 3 (yellow and blue colors) propagate with a speed equal to the sound speed of neutrals ($c_{\mathrm{n}0}/c_{\mathrm{tot}}=\sqrt{2}/2\approx 0.707$) for small values of α/k. For large α/k, these solutions have zero propagation speed (zero ω_R). The ratio ω_I/k is again small for weak or strong collisional coupling (small or large α/k), and has the maximum located on the x axis at a point very close to F_M.

These results are easy to interpret. For small collisional frequencies compared to the real part of the wave frequency, i.e., weak collisional coupling, the propagation of the waves mostly depends on the properties of either of the fluids. Since the neutrals and the charges have different sound speeds, in the case of weak collisional coupling the waves propagate either at the sound speed of neutrals or of charges. Correspondingly, if one perturbs charges with a certain amplitude, the weak drag forces will translate some of this perturbation to neutrals, but the amplitude of this perturbation will be very small, as indeed follows from the polarization relations in Eq. (30). The opposite is also true.

For large collision frequencies compared to the real part of the wave frequency, that is, strong collisional coupling, the charges and the neutrals become coupled and the wave propagates at the sound speed of the whole fluid. Subsequently, the amplitudes of the velocities of neutrals and charges are equal.

For the numerical tests, we only show the results corresponding to solution 4 (the green lines in Fig. 1), that is, the solution for waves which propagate at the sound speed of the charges for small α/k. The results in other cases are similar. We have selected one wavenumber so that αρ_tot/k/c_tot is less than F_M, and three values of k, for which this quantity is larger than F_M. These frequencies are marked with symbols in Fig. 1. In order to relate the results to observations, Table 2 provides the period, P = 2π/ω_R, and the damping time, T_D = 2π/ω_I for the waves used in the simulations (the wavenumbers k marked in Fig. 1), and for the wave corresponding to the maximum damping relative to the wavenumber (with wavenumber k_M) in physical units. We can see that the waves have extremely short periods and wavelengths and do not correspond to the waves observed in the chromosphere, which have typical frequencies of 3–5 mHz. This is due primarily to our use of a uniform, unstratified background atmosphere which does not correspond to the chromosphere. The temperatures of charges and neutrals differ by a factor of four, and this situation is also unrealistic. For these reasons, the value of the collisional parameter α corresponding to the background atmosphere may be unphysical for the chromosphere. Moreover, for the purposes of code verification we wanted to test both uncoupled and strongly coupled cases, and the values of k that we have chosen cover the full range, even if the waves corresponding to the largest values of k have small temporal and spatial scales that cannot be observed.

3.6. Solution of the dispersion relation for Alfvén waves

As in the case of the acoustic waves, we solve the dispersion relation for values of k between 6180 and 6.18 × 10⁻⁶ m⁻¹, and we plot the real and imaginary parts of the wave frequency in Fig. 2, as a function of the wavenumber in nondimensional units, C = f(D), defined in Eq. (40), i.e., ω_{I, R}/k/v_A0 as function of αρ_tot/k/v_A0. The dispersion relation is a third-order equation in ω and there are three solutions. In our case, ρ_c0/ρ_tot = 1/3, and the solutions obtained using Mathematica software are

$$\begin{array}{cc}& C_1=\frac{i}{3}(D+\frac{D^2-3}{G}+G),\hfill \\ \hfill & C_{3,2}=\frac{i}{3}(D+L_{3,2}^0{i}^{-4/3}\frac{D^2-3}{G}-L_{3,2}^1{i}^{-2/3}G),\hfill \end{array}$$(45)

with

$$\begin{array}{cc}& G={(D^3+3\sqrt{3-3D^2+D^4})}^{1/3},\hfill \\ \hfill & L_3^0=-i\sqrt{3}/2-1/2;L_2^0=1,\hfill \\ \hfill & L_3^1=i\sqrt{3}/2-1/2;L_2^1=1.\hfill \end{array}$$(46)

Fig. 2. Solutions of the dispersion relation, Eq. (39), for the Alfvén wave case. Left and right panels present the real and imaginary parts of the wave frequency ω, respectively, as a function of the wavenumber k. The axes are scaled in the nondimensional units: C and D, defined in Eq. (40) for better visualization. Three solutions are marked with different colors and line styles: red solid for solution 1, orange dashed for solution 2, and green dotted for solution 3. The values of k for which we compare the analytical and numerical solutions using the values of ω, corresponding to solution 3 (green dots), are marked in the panels: pink triangle (k = 61.8 m−1), green triangle (k = 6.18 m−1), orange star (k = 6.18 × 10−2 m−1), and black triangle (k = 6.18 × 10−6 m−1).

These represent one wave which does not propagate (ω_R = 0, red color), one wave traveling in the positive direction of z axis (green color), and a similar wave traveling in the negative direction (yellow color).

For small values of the collisional frequency compared to the real part of the wave frequency, the propagation speed (v_ph = ω_R/k) is the Alfvén speed considering only the density of charges (v_A0). For the large values of α/k, the value of v_ph is the Alfvén speed of the whole plasma considering both neutrals and charges, v_Atot. We note that ${{\mathit{v}}_{\mathrm{A}}}_{\mathrm{tot}}/{{\mathit{v}}_{\mathrm{A}}}_0=\sqrt{{\rho }_{\mathrm{c}0}/{\rho }_{\mathrm{tot}}}=1/\sqrt{3}\approx 0.577$. With no collisions between charges and neutrals, the imaginary part of solutions 2 and 3 vanishes,

$$\begin{array}{c}\hfill C_{3,2}=\pm 1,\end{array}$$(47)

meaning that the damping is related to the presence of neutrals. In our case, the imaginary part is positive, and the maximum damping relative to the wavenumber, the ratio ω_I/k, is located at a point $D_{\mathrm{M}}=\alpha {\rho }_{\mathrm{tot}}/k_{\mathrm{M}}^{\mathrm{A}}/{\mathit{v}}_{\mathrm{A}0}=2/\sqrt{3}$ on the x axis, corresponding to a value of $k_{\mathrm{M}}^{\mathrm{A}}\approx 8.8$ m⁻¹. This gives the value of ω_R(D = D_M)/(αρ_n0) ≈ 1. As in the case of the acoustic waves, the real part of the wave frequency is almost equal to the ion-neutral collision frequency at the point where the damping relative to wavenumber is maximum. Otherwise, the damping relative to the wavenumber is very small.

For the numerical tests below we use the solution traveling in the positive direction of the z axis (green color in Fig. 2) and for values of the wavenumbers, k, indicated in the legend of Fig. 2. The period and the damping time for the waves used in the simulations (corresponding to the values of k marked in Fig. 2), and for the wave corresponding to the maximum damping relative to the wavenumber (with wavenumber $k_{\mathrm{M}}^{\mathrm{A}}$) are given in Table 3. For the same reasons as the acoustic wave tests, the frequencies and the wavenumbers of the Alfvén waves do not correspond to typically observed waves in the chromosphere.

4. Results of numerical calculations

We ran simulations of acoustic and Alfvén waves using the equilibrium atmosphere and perturbation described in Sects. 3.3 and 3.4 as initial conditions. We used several values of the wavenumber k marked in the corresponding panels described in the above section. In order to observe the damping in time we choose a point located in the center of the domain and plot the evolution of the variables at the same point as function of time. The numerical solution is then compared to the analytical solution.

4.1. Temporal damping of acoustic waves

Figure 3 shows the numerical and analytical solutions for the acoustic waves with wavenumber values of k = 618, 6.18, 2.06, and 6.18 × 10⁻⁶ m⁻¹.

Fig. 3. Simulations of acoustic wave propagation in a homogeneous plasma. Numerically (red dashed line) and analytically (green solid line) calculated time evolution of the velocity of neutrals and charges as a function of time. Below the panels showing the velocity of neutrals, the difference (Δ) between the numerical solution and the analytical solution for unz is given. Time is measured in units of the wave period, 2π/ωR. Panels from left to right and from top to bottom: simulations for different values of the wavenumber k, as indicated in the legend of Fig. 1.

We can observe that for different values of k the amplitude damping and propagation velocity vary in agreement with what is expected from Fig. 1. For the large value of k = 618 m⁻¹ (panel a of Fig. 3) there is a phase shift between the oscillating fluid velocities of neutrals and charges, and both are slightly damped. We highlight the fact that the amplitude of oscillations of velocity of charges is significantly larger than that of the neutrals. For the intermediate values of k = 6.18 and 2.06 m⁻¹ (panels b and c of Fig. 3), we observe significant damping of both velocity oscillations but similar amplitudes and a smaller phase shift. For the small k case, that is, k = 6.18 × 10⁻⁶ m⁻¹ (panel d), the velocities of charges and neutrals are the same, and the wave propagation velocity is the sound speed of the plasma as a whole. This is expected since charges and neutrals are collisionally coupled.

We observe from Fig. 3 that numerical and analytical solutions are in very good agreement; in all panels, we plotted the difference between the numerical and the analytical solution, labeled Δ right below the plots for u_cx. The error can be seen to accumulate over time, but this is nevertheless very small.

The values from Table 2 are consistent with the values deduced from Fig. 1. The amplitude of a wave will decrease to a fraction f = exp(−2πP/T_D) in a period. For example, in the case of the wave with wavenumber k = 6.18 m⁻¹, f ≈ exp(−9.6/7.37) ≈ 0.27, which is consistent with the value observed in panel b of Fig. 3.

4.2. Temporal damping of Alfvén waves

Figure 4 shows the numerical and analytical solutions for the Alfvén waves with the wavenumber k set to the following values, k = 61.8, 6.18, 6.18 × 10⁻², and 6.18 × 10⁻⁶ m⁻¹. These values of k are marked in Fig. 2.

Fig. 4. Simulations of Alfvén wave propagation in a homogeneous plasma. Here we show the numerically (red dashed line) and analytically (green solid line) calculated time evolution of the velocity of neutrals and charges, and the magnetic field perturbation, as a function of time. Below the panels showing the velocity of charges, the difference (Δ) between the numerical solution and the analytical solution for ucx is plotted. Time is measured in units of the wave period, 2π/ωR. Panels from left to right and from top to bottom: simulations for different values of the wavenumber k, as indicated in the legend of Fig. 2.

We observe that numerical and analytical solutions match exactly and the amplitude damping and wave propagation speed vary with the ratio α/k as expected from Fig. 2.

For small α/k, the wave propagates at the Alfvén speed of the charges and the velocity of neutrals is almost zero (panel a of Fig. 4). For large α/k, the wave propagates at the Alfvén speed of the whole fluid and the velocities of neutrals and charges are equal (panel d). Yet again, for intermediate values of α/k, the damping of all perturbations is observed (panels b and c).

For all the panels in Fig. 4, we show the difference between the numerical and the analytical solution of u_cx, labeled Δ, right below the plots for u_cx. The error, similarly to the case of the acoustic waves, accumulates over time and is small. For the Alfvén waves, the values from Table 3 are also consistent with the values obtained from the numerical simulations. For example, the wave with wavenumber k = 6.18 m⁻¹ should decrease the amplitude to a fraction f ≈ exp(−1.63/2.76) ≈ 0.55 in a period, and this value can also be deduced from panel b of Fig. 4.

4.3. Time convergence test

In order to quantify the degree of agreement between the analytical and numerical solutions, we performed a time-step convergence test. We study the normalized error, ε(f, g), between analytical and numerical solutions as a function of the normalized integration time step, i.e. the ratio between the time step of the simulations, Δt, and the time step imposed by the single-fluid CFL condition, Δt_CFL. The normalized error is defined as

$$\begin{array}{c}\hfill \epsilon (f,g)=\sqrt{\frac{\sum _i{(f(x_i)-g(x_i))}^2}{\sum _i{(g(x_i))}^2}}·\end{array}$$(48)

In this expression, f and g are numerical and analytical solutions, respectively, computed at grid points x_i with i = 1, …N. The analytical solution, g, is independent of Δt. The normalized error defined this way was computed independently for all perturbed variables of the simulations. These are u_cz, u_nz, ρ_c1, ρ_n1, p_c1, and p_n1 for the acoustic waves, and u_nx, u_cx, and B_1x for the Alfvén waves. The results for different variables came to be very similar, therefore we only show the error for the velocity of charges, u_cz (acoustic waves) and u_cx (Alfvén waves). Computations for different wavenumber k also lead to very similar results. In the tests below we used k = 6.18 m⁻¹ for acoustic waves and k = 6.18 × 10⁻² m⁻¹ for Alfvén waves.

Figure 5 shows the results of the convergence tests for both types of waves. The errors are global errors computed after running simulations of acoustic and Alfvén waves until t = t_F, and the total time t_F corresponds to 15.96 periods for the acoustic waves and to 12.55 periods for the Alfvén waves. The maximum time step, Δt, used for the convergence test is close to the time step imposed by the single-fluid CFL condition. The numerical solutions were computed using time steps with a ratio of Δt/Δt_CFL= 0.812, 0.406, 0.203, 0.102, 0.051, and 0.025 for the acoustic waves and Δt/Δt_CFL= 0.84, 0.42, 0.21, 0.105, 0.053, and 0.026 for the Alfvén waves. The convergence tests were carried out for three values of the Toth et al. (2012) scheme parameter β of 0.5, 0.6 and 1.

Fig. 5. Normalized errors, ε, defined by Eq. (48), in the velocity of charges as a function of the ratio between the integration time step Δt and ΔtCFL, accumulated after running the simulation until t = tF. Panel a: simulations of acoustic waves, and panel b: simulations of Alfvén waves. Blue triangles show results for β = 1 (defined in Eq. (21)); green circles are for β = 0.6, and red crosses are for β = 0.5. Solid lines show the results of quadratic polynomial fit ε = a(Δt/ΔtCFL)2 + b(Δt/ΔtCFL)+c, with coefficients given in Tables 4 and 5.

Figure 5 shows that, as expected, the error between the analytical and the numerical solutions increases with the increase of the integration time step. However, this increase is different for different values of the parameter β. This parameter, defined in Eq. (21) determines the fraction of the implicit contribution to the numerical solution. For β = 1 the scheme is expected to be first-order accurate in time, and for β = 0.5 it is expected to be second-order accurate. In order to test if this is the case in our implementation, we performed a polynomial fit to the error function using the expression ε = a(Δt/Δt_CFL)² + b(Δt/Δt_CFL)+c, with the constraint a  ≥  0. The values of the coefficients together with their standard deviation resulting from the fit are given in Table 4 for the acoustic waves and in Table 5 for the Alfvén waves. The results confirm that the polynomial fit to the numerical values is dominantly second order for β ≈ 0.5 and first order for β = 1. We also observe that while higher-order convergence of numerical schemes is generally desirable, the absolute value of the error for the β = 1 tests at a practical time step is of similar magnitude (acoustic waves) or smaller (Alfvén waves) than for the β ≈ 0.5 tests.

It must be noted that the numerical solution here uses Δt up to the value dictated by the explicit single-fluid CFL condition. Therefore, we conclude that our implementation allows us to efficiently overcome the small time-step limitations implied by the stiff collisional terms in the two-fluid model.

5. Waves in a gravitationally stratified atmosphere

In this section we test the capabilities of the code to model waves in a strongly gravitationally stratified solar chromosphere. We assume a model atmosphere with all hydrodynamical parameters and purely horizontal magnetic field, B_x0, stratified in the vertical direction, z. The temperature is considered uniform (height-independent), and different for charges and neutrals. If we neglect viscosity, and consider only elastic collisions, adiabatic equation of energy, and ideal Ohm’s law, the linearized equations can be written as follows.

$$\begin{array}{cc}& \frac{\partial {\rho }_{\mathrm{c}1}}{\partial t}+u_{\mathrm{cz}}\frac{\mathrm{d}{\rho }_{\mathrm{c}0}}{\mathrm{d}z}+{\rho }_{\mathrm{c}0}\frac{\partial u_{\mathrm{cz}}}{\partial z}=0,\hfill \\ \hfill & \frac{\partial {\rho }_{\mathrm{n}1}}{\partial t}+u_{\mathrm{nz}}\frac{\mathrm{d}{\rho }_{\mathrm{n}0}}{\mathrm{d}z}+{\rho }_{\mathrm{n}0}\frac{\partial u_{\mathrm{nz}}}{\partial z}=0,\hfill \\ \hfill & {\rho }_{\mathrm{c}0}\frac{\partial u_{\mathrm{cz}}}{\partial t}+{\rho }_{\mathrm{c}1}g+\frac{\partial p_{\mathrm{c}1}}{\partial z}+\frac{1}{{\mu }_0}(\frac{\partial B_{\mathrm{x}1}}{\partial z}{B}_{\mathrm{x}0}+\frac{\mathrm{d}{B}_{\mathrm{x}0}}{\mathrm{d}z}{B}_{\mathrm{x}1})\hfill \\ \hfill & =\alpha {\rho }_{\mathrm{n}0}{\rho }_{\mathrm{c}0}(u_{\mathrm{nz}}-u_{\mathrm{cz}}),\hfill \\ \hfill & {\rho }_{\mathrm{n}0}\frac{\partial u_{\mathrm{nz}}}{\partial t}+{\rho }_{\mathrm{n}1}g+\frac{\partial p_{\mathrm{n}1}}{\partial z}=\alpha {\rho }_{\mathrm{n}0}{\rho }_{\mathrm{c}0}(u_{\mathrm{cz}}-u_{\mathrm{nz}}),\hfill \\ \hfill & \frac{\partial p_{\mathrm{c}1}}{\partial t}-c_{\mathrm{c}0}^2\frac{\partial {\rho }_{\mathrm{c}1}}{\partial t}-c_{\mathrm{c}0}^2{u}_{\mathrm{cz}}\frac{\mathrm{d}{\rho }_{\mathrm{c}0}}{\mathrm{d}z}+u_{\mathrm{cz}}\frac{\mathrm{d}{p}_{\mathrm{c}0}}{\mathrm{d}z}=0,\hfill \\ \hfill & \frac{\partial p_{\mathrm{n}1}}{\partial t}-c_{\mathrm{n}0}^2\frac{\partial {\rho }_{\mathrm{n}1}}{\partial t}-c_{\mathrm{n}0}^2{u}_{\mathrm{nz}}\frac{\mathrm{d}{\rho }_{\mathrm{n}0}}{\mathrm{d}z}+u_{\mathrm{nz}}\frac{\mathrm{d}{p}_{\mathrm{n}0}}{\mathrm{d}z}=0,\hfill \\ \hfill & \frac{\partial B_{\mathrm{x}1}}{\partial t}+B_{\mathrm{x}0}\frac{\partial u_{\mathrm{cz}}}{\partial z}+u_{\mathrm{cz}}\frac{\mathrm{d}{B}_{\mathrm{x}0}}{\mathrm{d}z}=0.\hfill \end{array}$$(49)

We separate the time dependence, assumed to be of the form exp(iωt) with constant ω, and combine the system into the equation of the vertical velocity for the charges: u_cz(z, t) = ũ_cz(z)exp(iωt), obtaining a fourth-order ODE:

$$\begin{array}{cc}& \frac{{\mathrm{d}}^4{\stackrel{\sim }{u}}_{\mathrm{cz}}}{\mathrm{d}{z}^4}{a}_{\mathrm{c}}{a}_{\mathrm{n}}+\frac{{\mathrm{d}}^3{\stackrel{\sim }{u}}_{\mathrm{cz}}}{\mathrm{d}{z}^3}(a_{\mathrm{n}}{b}_{\mathrm{c}}+a_{\mathrm{c}}{b}_{\mathrm{n}})\hfill \\ \hfill & +\frac{{\mathrm{d}}^2{\stackrel{\sim }{u_{\mathrm{c}}}}_z}{\mathrm{d}{z}^2}[b_{\mathrm{c}}{b}_{\mathrm{n}}+{\omega }^2(a_{\mathrm{c}}+a_{\mathrm{n}})-i\alpha \omega (a_{\mathrm{c}}{\rho }_{\mathrm{c}0}+a_{\mathrm{n}}{\rho }_{\mathrm{n}0})]\hfill \\ \hfill & +\frac{\mathrm{d}{\stackrel{\sim }{u}}_{\mathrm{cz}}}{\mathrm{d}z}\omega [\omega (b_{\mathrm{c}}+b_{\mathrm{n}})-i\alpha (b_{\mathrm{c}}{\rho }_{\mathrm{c}0}+b_{\mathrm{n}}{\rho }_{\mathrm{n}0})]\hfill \\ \hfill & +u_{\mathrm{cz}}{\omega }^3[\omega -i\alpha ({\rho }_{\mathrm{c}0}+{\rho }_{\mathrm{n}0})]=0.\hfill \end{array}$$(50)

where

$$\begin{array}{cc}\hfill a_{\mathrm{c}}(z)& ={c_{\mathrm{c}0}}^2+{v_{A0}}^2,a_{\mathrm{n}}(z)=c_{\mathrm{n}0}^2,b_{\mathrm{c},\mathrm{n}}(z)=\frac{1}{{\rho }_{\mathrm{c},\mathrm{n}0}}\frac{\partial \left({\rho }_{\mathrm{c},\mathrm{n}0}{a}_{\mathrm{c},\mathrm{n}}\right)}{\partial z},\hfill \\ \hfill c_{\mathrm{c},\mathrm{n}0}^2& =\gamma \frac{p_{\mathrm{c},\mathrm{n}0}}{{\rho }_{\mathrm{c},\mathrm{n}0}},\text{and}{v_{A0}}^2=\frac{{B_{\mathrm{x}0}}^2}{{\mu }_0{\rho }_{\mathrm{c}0}}·\hfill \end{array}$$(51)

The equilibrium variables for neutrals and for charges must fulfill the equations of state: Eq. (13) and the hydrostatic (HS) and magneto-hydrostatic (MHS) equilibrium conditions, respectively:

$$\begin{array}{cc}& \frac{\mathrm{d}{p}_{\mathrm{n}0}}{\mathrm{d}z}+{\rho }_{\mathrm{n}0}g=0,\hfill \\ \hfill & \frac{\mathrm{d}}{\mathrm{d}z}(p_{\mathrm{c}0}+{B_{\mathrm{x}0}}^2/(2{\mu }_0))+{\rho }_{\mathrm{c}0}g=0.\hfill \end{array}$$

Since the temperature is assumed uniform, the pressure for neutrals has an exponential profile with a uniform scale height, and the sound speed of neutrals is constant. If we consider that the magnetic pressure has the same scale height as the pressure on the charges, after solving HS/MHS equations we obtain

$$\begin{array}{cc}\hfill p_{\mathrm{n}0}(z)& =p_{\mathrm{n}0}(z_0)exp(-\frac{z}{H_{\mathrm{n}}}),\hfill \\ \hfill p_{\mathrm{c}0}(z)& =p_{\mathrm{c}0}(z_0)exp(-\frac{z}{H_{\mathrm{c}}}),\hfill \\ \hfill B_{\mathrm{x}0}(z)& =B_{\mathrm{x}0}(z_0)exp(-\frac{z}{2H_{\mathrm{c}}}),\hfill \end{array}$$(52)

with uniform scale heights:

$$\begin{array}{cc}\hfill H_{\mathrm{n}}& =\frac{k_{\mathrm{B}}{T}_{\mathrm{n}0}}{m_{\mathrm{H}}g},\hfill \\ \hfill H_{\mathrm{c}}& =\frac{2k_{\mathrm{B}}{T}_{\mathrm{c}0}}{m_{\mathrm{H}}g}\frac{p_{\mathrm{c}0}(z_0)+B_{\mathrm{x}0}{(z_0)}^2/(2{\mu }_0)}{p_{\mathrm{c}0}(z_0)}·\hfill \end{array}$$(53)

The densities obtained from the ideal gas laws for neutrals and charges (Eq. (13)) also have an exponential profile. In these conditions, the quantities defined in Eq. (51), namely a_c, n, b_c, n = −a_c, n/H_c, nn and α, are uniform; however, there are coefficients that explicitly contain density. We assumed that the scale of the height variation of the nonuniform coefficients in Eq. (50) is large compared to the oscillation wavelength, and therefore we can search for the solution in terms of plane waves, as

$$\begin{array}{c}\hfill {\stackrel{\sim }{u}}_{\mathrm{cz}}(z)=V_{\mathrm{c}}\text{exp}(-ikz),\end{array}$$

with V_c being the uniform real amplitude, and k the uniform complex wavenumber. The other variables are also in the form of plane waves (including the time dependence):

$$\begin{array}{cc}\hfill \{u_{\mathrm{nz}},\frac{{\rho }_{\mathrm{c}1}}{{\rho }_{\mathrm{c}0}},\frac{{\rho }_{\mathrm{n}1}}{{\rho }_{\mathrm{n}0}},\frac{p_{\mathrm{c}1}}{p_{\mathrm{c}0}},\frac{p_{\mathrm{n}1}}{p_{\mathrm{n}0}},\frac{B_{\mathrm{x}1}}{B_{\mathrm{x}0}}\}& =\{{\stackrel{\sim }{V}}_{\mathrm{n}},{\stackrel{\sim }{R}}_{\mathrm{c}},{\stackrel{\sim }{R}}_{\mathrm{n}},{\stackrel{\sim }{P}}_{\mathrm{c}},{\stackrel{\sim }{P}}_{\mathrm{n}},\stackrel{\sim }{B}\}\hfill \\ \hfill & \times exp(i(\omega t-kz)),\hfill \end{array}$$(54)

where Ṽ_n, R̃_c, R̃_n, P̃_c,P̃_n, and B̃ are uniform complex amplitudes. The following relations are obtained,

$$\begin{array}{cc}& {\stackrel{\sim }{V}}_{\mathrm{n}}=V_{\mathrm{c}}(1+\frac{i({\omega }^2-ik{b}_{\mathrm{c}}-k^2{a}_{\mathrm{c}})}{\alpha \omega {\rho }_{\mathrm{n}0}}),\hfill \\ \hfill & {\stackrel{\sim }{R}}_{\mathrm{c},\mathrm{n}}=\frac{{\stackrel{\sim }{V}}_{\mathrm{c},\mathrm{n}}}{\omega }(k-i\frac{1}{H_{\mathrm{c},\mathrm{n}}}),\hfill \\ \hfill & {\stackrel{\sim }{P}}_{\mathrm{c},\mathrm{n}}=\frac{{\stackrel{\sim }{V}}_{\mathrm{c},\mathrm{n}}}{\omega }(k\gamma -i\frac{1}{H_{\mathrm{c},\mathrm{n}}}),\hfill \\ \hfill & \stackrel{\sim }{B}=\frac{V_{\mathrm{c}}}{\omega }(k-i\frac{1}{2H_{\mathrm{c}}}).\hfill \end{array}$$(55)

The dispersion relation is a fourth-order equation in ω:

$$\begin{array}{cc}& (-{\omega }^2+k^2{a}_{\mathrm{c}}+ik{b}_{\mathrm{c}}+i\omega \alpha {\rho }_{\mathrm{n}0})\hfill \\ \hfill & \times (-{\omega }^2+k^2{a}_{\mathrm{n}}+ik{b}_{\mathrm{n}}+i\omega \alpha {\rho }_{\mathrm{c}0})+{\omega }^2{\alpha }^2{\rho }_{\mathrm{n}0}{\rho }_{\mathrm{c}0}=0.\hfill \end{array}$$(56)

In our particular case, the wavenumber k obtained from this space-dependent dispersion relation turns out to be almost uniform, and therefore we consider the plane wave solution to be a fair approximation.

5.1. Initial conditions

5.1.1. Equilibrium atmosphere

We choose the temperature for neutrals to be T_n0 = 6000 K. Further, for the number density of the neutrals and charges we use the values taken from the VALC (Vernazza et al. 1981) model at z₀ ≈ 500 km: n_c0(z₀) = 5 × 10¹⁷ m⁻³, and n_n0(z₀) = 2.1 × 10²¹ m⁻³. In this test we take the temperature of the electrons into account; unlike the tests in the uniform atmosphere, they contribute to the pressure of charges and to collisions. We choose the value of B_x0(z₀) = 10⁻⁴ T. In order to have the same scale height for neutrals and charges, we use the following temperature for the charges.

$$\begin{array}{c}\hfill T_{\mathrm{c}0}=\frac{n_{\mathrm{c}0}(z_0)k_{\mathrm{B}}{T}_{\mathrm{n}0}-B_{\mathrm{x}0}{(z_0)}^2/(2{\mu }_0)}{2n_{\mathrm{c}0}(z_0)k_{\mathrm{B}}}.\end{array}$$(57)

This gives the value of T_c0 ≈ 2422.47 K and H_n = H_c ≈ 1.8 × 10⁵ m. Subsequently, the pressures of charges and neutrals at z₀, p_c0(z₀) and p_n0(z₀) are obtained from the ideal gas law (Eq. (13)). In these conditions, the wavenumber k obtained from the dispersion relation Eq. (56) is almost uniform, and has a value of k ≈ 1.383 × 10⁻⁴ + 2.767 × 10⁻⁶i m⁻¹. Therefore, the wavelength is about four times shorter than the density scale height. The densities are calculated afterwards from the ideal gas laws for neutrals and charges (Eq. (13)), taking ρ_n0 = n_nm_H, ρ_c0 = n_em_H. We cover the domain L_z = 1.6 Mm with 32 000 grid points.

5.1.2. Perturbation

We choose the period of the wave to be P = 5 s, and calculate the frequency (ω = 2π/P) and the wavenumber (k) from the dispersion relation, Eq. (56). We choose the amplitude of the perturbation of the velocity of charges as a fraction of the background sound speed, V_c = 10⁻³c₀, where

$$\begin{array}{c}\hfill c_0=\sqrt{\gamma (p_{\mathrm{n}0}+p_{\mathrm{c}0})/({\rho }_{\mathrm{n}0}+{\rho }_{\mathrm{c}0})}\end{array}$$(58)

is the sound speed of the whole fluid, and its value is c₀ ≈ 9.1 km s⁻¹. We calculate the amplitudes of the other perturbations from the polarization relations, Eq. (55). The perturbation is generated by a driver at each time step in the ghosts points at the bottom of the atmosphere, which makes it the lower boundary condition. The perfectly matched layer (PML) is used as the upper boundary condition and is specially designed to absorb waves without reflections. It was first introduced for the first time for electromagnetic waves in the Maxwell equation by Berenger (1994), applied to Euler equations by Hu (1996), and to acoustic waves in a strongly stratified solar convection zone by Parchevsky & Kosovichev (2007). The description of the implementation of the PML in the Mancha code can be found in Felipe et al. (2010). In this test we have used the value of the scheme parameter β = 1.

5.2. Results

We run the code in two regimes. In the first case, we solved fully nonlinear equations for perturbations, and in the second case we evolved linearized equations where only the first-order terms were kept. The latter was done for comparison purposes, since the analytical solution assumes linear regime. We compare the numerical solutions at time t = 215.115 s with the analytical solution. At this time, simulations reached the stationary state, since the wave has reached the upper boundary and several periods of the wave have passed through the boundary. The results are shown in Figs. 6 and 7.

Fig. 6. Analytical solution (green solid line) and numerical linear solution (red dashed line) for the vertical velocity of charges (uzc) and neutrals (uzn), and for the perturbation in the x component of the magnetic field (B1x) for a snapshot taken in stationary state at t = 215.115 s. Below the panel showing the solutions for uzc, the difference between the numerical linear solution and the analytical solution for uzc is plotted.

Fig. 7. First panel: analytical solution (green solid line) and numerical solutions: linear (red dashed line) and nonlinear (blue dotted line) for uzc for a snapshot taken in stationary state at t = 215.115 s. Second and third panels: decoupling (uzn − uzc) for the analytical and linear numerical solutions, and for the nonlinear numerical solutions, respectively.

Figure 6 shows the analytical solution (green, solid line) superposed on the linear numerical solution (red, dashed line) for the vertical velocity of charges (u_zc) in the first panel, the vertical velocity of neutrals (u_zn) in the second panel, and for the perturbation in the x component of the magnetic field (B_x1) in the third panel. Below the panel of u_zc we show the difference (Δ) between the linear numerical and the analytical solution of u_zc. We observe that the numerical linear solution is in very good agreement with the analytical solution for the three quantities considered, and that the error is small (below 2%).

The first panel of Fig. 7 shows the nonlinear effects by superposing the nonlinear numerical solution (blue, dotted line) of u_zc along with the analytical solution (green solid line) and the linear numerical solution (red dashed line), for the same snapshot taken at t = 215.115 s. In the second and third panels we show the decoupling in the vertical velocity (u_nz − u_cz) for the analytical and linear solutions superposed, and the nonlinear solution, respectively.

We observe that the wave profile steepens at the end of the domain when the amplitude becomes large and nonlinear effects are visible in the case of the nonlinear solution. As a consequence, the amplitude of the wave is smaller than in the linear case (see e.g., Landau & Lifshitz 1987). Wave amplitude grows with height in a gravitationally stratified atmosphere because of the density decrease. In the case considered here, the damping is not large enough to overcome this growth,and therefore the wave evolves into a shock. The amplitude growth is nevertheless below the growth in the MHD limit.

We also observe that the decoupling predicted analytically from the relation between Ṽ_n and V_c in Eq. (55) agrees with the decoupling obtained from the linear numerical simulation. In order to intuitively understand the reason for the linear decoupling, we show in Fig. 8 the relevant frequencies corresponding to this problem. Even though the collisional parameter α is uniform because of the uniform background temperature of neutrals and charges, and has the value ≈4.1 × 10¹¹ m³ kg⁻¹ s⁻¹, the collision frequency also depends on density (Eq. (5)), and has an exponential profile. While for the charges the collision frequency ν_cn is larger than both the ion-cyclotron frequency (ω_ci) and the wave frequency (which determine the hydrodynamical timescale), for the neutrals the neutral-charges collision frequency goes from being greater than the wave frequency to being less than the wave frequency at a point located at z ≈ 1.4 Mm in the atmosphere. This is the point after which we observe the decoupling for the linear numerical solution and the analytical solution. The code efficiently captures the transition from a coupled to a partially decoupled regime. In the nonlinear case, the hydrodynamical timescale is smaller than in the linear case. The shock front width, which determines the hydrodynamical space scale in the nonlinear case, is smaller than the wavelength. The nonlinear decoupling appears spatially at the shock front, and is almost five times larger than the linear decoupling.