© ESO 2019

1. Introduction

The solar chromosphere is a very dynamic region that is constantly perturbed by propagating waves and shocks. Waves can be observed by measuring fluctuations in the radiation intensity due to fluctuations in the density and temperature of both compressible acoustic and magneto-acoustic waves, and by measuring the Doppler shift associated with fluctuations in velocity for incompressible waves. Waves have been observed at different layers in the solar atmosphere, from the photosphere to the corona, but limited spatial and temporal resolution restricts their detection to low-frequency waves. Significant effort has been made to measuring waves in the solar corona (Verwichte et al. 2004; Ofman et al. 2008; Gupta 2014). Photospheric and chromospheric oscillations in different magnetic structures have been measured for decades. Making interpretations about chromospheric measurements is especially difficult because of the fast evolving structures, supersonic motions, and high-frequency oscillatory phenomena. Physical processes in the chromosphere happen at characteristic times dictated by the Alfvén speed that can reach several hundreds of km s⁻¹. A review of multi-wavelength studies of waves in the chromosphere is provided by Jess et al. (2015).

The following two main mechanisms are thought to be responsible for the heating of the solar atmosphere: magnetic reconnection and waves. These two mechanisms are not independent since the magnetic reconnection can be a source of waves (Verwichte et al. 2004; Jess et al. 2008; Luna et al. 2008; Li & Zhang 2012; Provornikova et al. 2018) and since waves can produce instabilities needed to trigger reconnection processes (Isobe & Tripathi 2006; Lee et al. 2014; McLaughlin et al. 2009). Even though the chromosphere has a much lower temperature than the solar corona, it has much higher energy density and requires more energy to compensate for the radiative losses (Jess et al. 2015; Withbroe & Noyes 1977; Anderson & Athay 1989). Therefore, mechanisms providing energy to the chromosphere have to be very efficient.

Contrary to previous studies, based on low resolution observations or on 1D simulations, high-frequency waves may play an important role in the energization of the solar chromosphere. From an observational point of view, waves with frequencies above 10–50 mHz are extremely difficult to detect and are often referred to as high-frequency waves (e.g., see Fossum & Carlsson 2006; Bello González et al. 2009). Musielak et al. (1994) considered wave generation in the solar convection zone and computed the acoustic wave energy flux as a function of frequency. Their theoretical calculations revealed that acoustic flux has a broad maximum in the high-frequency region around 100 mHz (periods ≈10 s). However, firm detection of the high-frequency part of the solar oscillation spectrum is still uncertain. For example, Bello González et al. (2009, 2010a,b), using high resolution data obtained by Gregor Fabry-Pérot Interferometer (GFPI) at Vacuum Tower Telescope (VTT) and Imaging Magnetograph eXperiment (IMaX) instrument aboard the Sunrise mission, revealed a high acoustic flux in the chromosphere for high-frequency waves. This result is in contradiction to the earlier work by Fossum & Carlsson (2006) that uses Transition Region and Coronal Explorer (TRACE) data.

The maxima of the observed solar wave spectra fall into 3–5 mHz range, depending on the height. This low-frequency part of the solar wave spectrum has been extensively studied both theoretically and observationally. From a theoretical point of view, a magnetohydrodynamic (MHD) approximation can be safely used for the studies of low-frequency wave phenomena. A review of theoretical studies of wave propagation through photospheric, chromospheric, and coronal structures can be found in Khomenko & Calvo Santamaria (2013), Khomenko & Collados (2015). Waves are sensitive to the horizontal and vertical structure of the atmosphere. In particular, strong vertical gradients in the acoustic and Alfvén speed lead to phenomena of mode transformation (Cally 2006; Cally & Goossens 2008) or fast magneto-acoustic mode refraction (Bogdan et al. 2003; Khomenko & Collados 2006) with particularly interesting phenomena happening around singular magnetic field structures, such as null points (Lee et al. 2014; Santamaria et al. 2017; Provornikova et al. 2018).

When considering high-frequency waves and shocks, damping or dissipation effects become important. Small-scale disturbances can easily be damped through various non-ideal mechanisms such as viscosity, radiation, thermal conduction, electrical resistivity, or ion-neutral friction. Other mechanisms related to the atmospheric structure are also at play, such as resonant absorption or phase mixing. The importance of these effects depends on the height in the solar atmosphere. For example, damping due to thermal conduction has been studied under coronal conditions (Klimchuk et al. 2004; Gupta 2014; Soler & Terradas 2015), phase mixing, and resonant absorption and are thought to play a role in damping coronal loop oscillations (Vasquez 2005; Verwichte et al. 2004). High-frequency compressible waves are heavily damped by radiative processes at the base of the photosphere where the energy that is dissipated by shocks is converted into radiation (Carlsson & Stein 2002). Viscous damping of shock waves is reconsidered by Arber et al. (2016). Additionally, resistive dissipation of Alfvén waves, enhanced through ion-neutral interaction, is considered by Song & Vasyliūnas (2011), Shelyag et al. (2016), Martínez-Sykora et al. (2016).

Under the conditions of the solar chromosphere, the reduction in the collisional coupling with height leads to partial decoupling of the neutral and charged components of the solar plasma, and produces a series of multi-fluid effects on waves (Ballester et al. 2018a). In order to treat these effects, classical MHD modeling may be insufficient and extensions to this approach are needed, such as the use of single-fluid formalism together with the generalized Ohm’s law, or a two-fluid (or multi-fluid) formalism. Ion-neutral effects on waves are extensively studied analytically under the assumption of homogeneous, unbounded plasma (Zaqarashvili et al. 2011, 2013; Soler et al. 2013a,b; Martínez-Gómez et al. 2016, 2017; Ballester et al. 2018b), or by using simplified magnetic configurations, such as flux tubes (Zaqarashvili et al. 2013; Soler et al. 2017). The main conclusions from these works is that ion-neutral effects produce important wave damping when the wave frequency is close to the ion-neutral collision frequency. In addition, new wave modes, that do not exist in a pure proton-electron plasma, can appear while cut-off wave numbers may appear for other modes.

Complex nonlinear dynamic processes and the 3D structure of the chromosphere strongly limit the practicality of linear theory for interpretation of observations, and numerical modeling techniques must be applied. Such numerical models mainly use a single-fluid approach that introduce the partial ionization effects through a generalized Ohm’s law (Khomenko & Collados 2012; Cheung & Cameron 2012; Martínez-Sykora et al. 2012, 2016; Khomenko et al. 2014a, 2017, 2018; Shelyag et al. 2016). One of the advantages of the numerical approach is the possibility to study wave dissipation and heating due to ion-neutral interactions. The heating is expressed via the ambipolar term in the single-fluid approach (Song & Vasyliūnas 2011; Khomenko 2017), or via the frictional heating term in the two-fluid approach (Leake et al. 2014), and is dropped in the linear theory, being quadratic in the perturbed quantities. Nonlinear simulations allow us to fully take the effects of wave heating due to neutrals into account. Most of such modeling has been done using the single-fluid approach, which is most appropriate for the dense layers of the solar atmosphere (photosphere and low chromosphere). Only recently has modeling been performed in the solar atmosphere for waves in the two-fluid fully nonlinear assumption (Maneva et al. 2017).

One particular aspect of the chromospheric physics that has been scarcely studied is the influence of multi-fluid effects on the formation and dissipation of chromospheric shock waves. Hillier et al. (2016) study the formation and evolution of slow-mode shocks driven by reconnection in a partially ionized plasma. A complex multi-fluid structure of the shock transition was modeled revealing structures similar to C-shocks or J-shocks in the classification by Draine & McKee (1993). The frictional heating associated with ion-neutral decoupling at the shock front was up to 2% of the available magnetic energy.

The aim of the present paper is to deepen our understanding of the multi-fluid nature of chromospheric shocks and their contribution to the chromospheric heating. To do so, we performed two-fluid modeling of propagation of fast magneto-acoustic waves generated at the base of the chromosphere for a particular case with a horizontal magnetic field. Depending on the conditions, these waves steepen into shocks at the middle-to-upper chromosphere, leading to decoupling between ion and neutral fluid velocities. Unlike many studies mentioned above, we fully take nonlinear effects into account. The simulations are performed using a newly extended MANCHA3D code (Popescu Braileanu et al. 2019). We study the effects of the wave frequency and amplitude, as well as the background magnetic field strength on the height where ions and neutrals become decoupled. The results of the numerical simulations are compared to analytical solutions within the two-fluid and single-fluid descriptions, thus allowing us to separate the nonlinear, linear, and two-fluid effects from wave damping and dissipation.

2. Numerical model

In the present study we used a two-fluid model assuming a purely hydrogen plasma. This allowed us to greatly simplify the expressions for the interaction terms in the two-fluid equations (Leake et al. 2014; Khomenko et al. 2014b). Nevertheless, it has to be noted that modeling based on hydrogen plasma cannot be consistently applied through atmospheric layers. The solar atmosphere is composed of multiple chemical species in different ionization states. Thermodynamic conditions in the photosphere are such that the main donors of electrons are metals, and therefore hydrogen plasma modeling cannot self-consistently take this feature into account. For this reason, modeling performed in this work only applies to the atmospheric layers starting from the upper photosphere and upward. Below, we describe the equations solved by the code, the equilibrium atmospheric model, and the wave driver.

2.1. Two-fluid equations

The following set of equations is solved numerically by the two-fluid version of the code MANCHA3D (Popescu Braileanu et al. 2019). Electrons and ions are evolved together by a single set of equations for charges (marked by subindex “c”) while neutrals are evolved separately (marked by subindex “n”). Here we neglected Hall, battery, and other smaller effects from the generalized Ohm’s law written in the system of reference of charges. Also, we did not consider ionization and recombination processes, viscosity, thermal conduction, or radiation. The only non-ideal terms are neutral-charges, elastic collisional terms in the momentum, and energy equations,

$$\begin{array}{cc}& \frac{\partial {\rho }_{\mathrm{n}}}{\partial t}+\mathrm{\nabla }·({\rho }_{\mathrm{n}}{\mathit{u}}_{\mathrm{n}})=0,\hfill \\ \hfill & \frac{\partial {\rho }_{\mathrm{c}}}{\partial t}+\mathrm{\nabla }·({\rho }_{\mathrm{c}}{\mathit{u}}_{\mathrm{c}})=0,\hfill \\ \hfill & \frac{\partial ({\rho }_{\mathrm{n}}{\mathit{u}}_{\mathrm{n}})}{\partial t}+\mathrm{\nabla }·({\rho }_{\mathrm{n}}{\mathit{u}}_{\mathrm{n}}{\mathit{u}}_{\mathrm{n}}+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}}{\mathit{u}}_{\mathrm{c}}+p_{\mathrm{c}})=\mathit{J}\times \mathit{B}+{\rho }_{\mathrm{c}}\mathit{g}-{\mathit{R}}_{\mathrm{n}},\hfill \\ \hfill & \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)+p_{\mathrm{n}}{\mathit{u}}_{\mathrm{n}})={\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)+p_{\mathrm{c}}{\mathit{u}}_{\mathrm{c}})={\mathit{u}}_{\mathrm{c}}·(\mathit{J}\times \mathit{B})\hfill \\ \hfill & +{\rho }_{\mathrm{c}}{\mathit{u}}_{\mathrm{c}}·\mathit{g}-M_{\mathrm{n}},\hfill \\ \hfill & \frac{\partial \mathit{B}}{\partial t}-\mathbf{\nabla }\times ({\mathit{u}}_{\mathrm{c}}\times \mathit{B})=0,\hfill \end{array}$$(1)

with

$$\begin{array}{c}\hfill {\mathit{R}}_{\mathrm{n}}=\alpha {\rho }_{\mathrm{n}}{\rho }_{\mathrm{c}}({\mathit{u}}_{\mathrm{c}}-{\mathit{u}}_{\mathrm{n}}),\end{array}$$(2)

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

In these equations, u_c, n, ρ_c, n, p_c, n, e_c, n, and T_c, n are center of mass velocities, mass densities, pressures, internal energies, and temperatures of charges and neutrals, B is the magnetic field, J is the current density, and R_n and M_n are the neutral momentum and energy collisional source terms, respectively. We neglected the contribution of electrons to the mass density of charges with m_e/m_H <  10⁻³, ρ_c = m_Hn_e + m_en_e ≈ m_Hn_e, with n_e being the electron number density. The gravitational acceleration g is oriented in the negative z direction. The collisional parameter α is defined as follows,

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

where ν_en and ν_in are the electron-neutral and ion-neutral collision frequencies, respectively,

$$\begin{array}{c}\hfill {\nu }_{\mathrm{in}}=n_{\mathrm{n}}\sqrt{\frac{16k_{\mathrm{B}}{T}_{\mathrm{cn}}}{\pi m_{\mathrm{H}}}}{\sigma }_{\mathrm{in}};{\nu }_{\mathrm{en}}=n_{\mathrm{e}}\sqrt{\frac{8k_{\mathrm{B}}{T}_{\mathrm{cn}}}{\pi m_{\mathrm{e}}}}{\sigma }_{\mathrm{en}};\end{array}$$(5)

with T_cn being the average temperature of neutrals and charges, σ_in = 5 × 10⁻¹⁹ m² and σ_en = 10⁻¹⁹ m² (Braginskii 1965).

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

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

The factor of two in the last equation appears to be due to the contributions of both ions and electrons. Additionally, both temperatures are assumed to be equal to T_c.

The equations were advanced in time using a semi-implicit scheme from Tóth et al. (2012) with the collisional terms implemented implicitly, and the rest of the terms were implemented in a second order accurate Runge Kutta scheme. The details of the implementation are provided in Popescu Braileanu et al. (2019). The collisional timescales can be significantly smaller than the hydrodynamical timescale, making the numerical discretization of the equations stiff. The implicit implementation of the collisional terms is not affected by a very small time step limitation imposed by the collisional time scale. This makes it possible to use the time step determined by the ideal MHD CFL condition and thus advance in time about 75 times faster as compared to a fully explicit implementation.

2.2. Equilibrium atmosphere

We assumed a model atmosphere with all hydrodynamic variables and a purely horizontal magnetic field, B_x0, that is stratified in the vertical, z, direction, and subindex “0” refers to equilibrium variables. The MANCHA3D code requires variables for the equilibrium atmosphere to be set separately for the neutral and charged components. In practice, any set of equilibrium variables can be used as long as they fulfill the conditions of purely hydrostatic (HS) equilibrium for neutrals, and magneto-hydrostatic (MHS) equilibrium for charges.

In equilibrium, we consider that the charges and neutrals have the same background temperature and we used the temperature profile from the model proposed by Vernazza, Avrett, and Loeser for the photosphere and Chromosphere (VALC; Vernazza et al. 1981), starting from heights, z, above z₀ ≈ 500 km, and ending just below the transition region at a temperature of 9000 K, at z_f ≈ 2.1 Mm. Apart from the temperature structure, no other stratification from the VALC model can be used directly because the equations mentioned above are only valid for hydrogenate plasma. Therefore, we recalculated the stratification of pressure and density of both neutrals and charges, assuming pure hydrogen, in order to satisfy the requirements of the HS and MHS equilibrium as explained above.

Apart from temperature, the only other parameter taken from VALC is the number density of charges and neutrals at the reference base level, z₀ ≈ 500 km, with n_c0(z₀) = 5 × 10¹⁷ m⁻³, and n_n0(z₀) = 2.1 × 10²¹ m⁻³. Then, the pressures of charges and neutrals at z₀, p_c0(z₀) and p_n0(z₀) are obtained from the ideal gas law from Eq. (6). As discussed above, we only considered the height range z = {z₀, z_f} because the hydrogen plasma model does not provide correct electron number densities in the deeper layers as most available electrons come from metals.

Taking Eq. (6) into account, the equations of HS and MHS balance for neutrals and charges as follows,

$$\begin{array}{c}\hfill \frac{\mathrm{d}{p}_{\mathrm{n}0}}{\mathrm{d}z}=-\frac{m_{\mathrm{H}}g}{k_{\mathrm{B}}T}{p}_{\mathrm{n}0},\end{array}$$(7)

$$\begin{array}{c}\hfill \frac{\mathrm{d}(p_{\mathrm{c}0}+p_{\mathrm{m}0})}{\mathrm{d}z}=-\frac{m_{\mathrm{H}}g}{2k_{\mathrm{B}}T}{p}_{\mathrm{c}0},\end{array}$$(8)

where p_m0 is the magnetic pressure, which is the only contribution from the Lorentz force produced by the horizontal magnetic field. The HS equation for neutrals is readily integrated, leading to:

$$\begin{array}{c}\hfill p_{\mathrm{n}0}(z)=p_{\mathrm{n}0}(z_0)exp(-\frac{m_{\mathrm{H}}g}{k_{\mathrm{B}}}{{\displaystyle \int }}_{z_0}^z\frac{1}{T(z^{\prime })}\mathrm{d}{z}^{\prime }).\end{array}$$(9)

In order to be able to integrate the MHS equation for charges, we assumed that the charges pressure and the magnetic pressure have the same exponential dependence on height,

$$\begin{array}{cc}& p_{\mathrm{c}0}(z)=p_{\mathrm{c}0}(z_0)exp(-2F(z))\hfill \end{array}$$(10)

$$\begin{array}{cc}& p_{\mathrm{m}0}(z)=\frac{B_{x0}{(z)}^2}{2{\mu }_0}=\frac{B_{x0}{(z_0)}^2}{2{\mu }_0}exp(-2F(z))+C\hfill \end{array}$$(11)

where F should satisfy F(z) ≥ 0 and F(z₀) = 0, and C is a constant. In introducing these expressions into Eq. (8), the following expression for F(z) was obtained,

$$\begin{array}{c}\hfill F(z)=\frac{m_{\mathrm{H}}g}{4k_{\mathrm{B}}}\frac{p_{\mathrm{c}0}(z_0)}{p_{\mathrm{c}0}(z_0)+B_{x0}{(z_0)}^2/2{\mu }_0}{{\displaystyle \int }}_{z_0}^z\frac{1}{T(z^{\prime })}\mathrm{d}{z}^{\prime }.\end{array}$$(12)

By setting the value of B_x0(z₀) = 10⁻⁴ T, a pressure profile of the charged component similar to that of the VALC model was obtained. Given the temperature and pressure profiles with height, the densities were calculated afterward from the ideal gas laws for neutrals and charges, Eqs. (6), taking ρ_n0 = n_nm_H, ρ_c0 = n_em_H.

The integration constant C in Eq. (11) allowed us flexibility in selecting the strength of the resulting magnetic field. We set,

$$\begin{array}{c}\hfill C=p_{\mathrm{n}0}{|}_{z=z_0+j{L}_z}-\frac{B_{x0}{(z_0)}^2}{2{\mu }_0}exp(-2F{|}_{z=z_0+j{L}_z}),\end{array}$$(13)

and B_x0(z) was then recovered from Eq. (11). In this equation, j ∈ [0, 1] is a fraction of the total vertical domain length, L_z = 1.6 Mm. By selecting j inside the domain, we made sure that neutral pressure and magnetic pressure were equal to each other at z = jL_z, thus creating two zones in the atmosphere with either neutral gas pressure or magnetic pressure that is more prevalent. In the following, we used two magnetic field profiles, S-profile with j = 1/2 and B-profile with j = 0. In the case of the S profile, the average field strength is smaller, about 1.5 × 10⁻³ T, while for the B-profile, it reaches 1.74 × 10⁻² T. We note that the integration constant dominates over the space varying term, having the value of 0.988 Pa for the S profile and 120.73 Pa for the B profile, thus making the magnetic field profiles almost flat. The difference between the maximum value of the magnetic field, obtained at the base of the atmosphere and the minimum value of the magnetic field obtained at the top of the atmosphere, is maximum for the minimum magnetic field profile that fulfills MHS. This happens when the integration constant is chosen as

$$\begin{array}{c}\hfill -\frac{B_{x0}{(z_0)}^2}{2{\mu }_0}exp(-2F{|}_{z=L_z}).\end{array}$$(14)

In this case, the magnetic field is zero at the top of the atmosphere and has a value of less than 1G at the base of the atmosphere.

The resulting stratification of pressures, number densities, temperature, and the magnetic field are displayed in Fig. 1. It can be observed in panel a that the neutral and magnetic pressures become equal in the middle of the domain for the S profile, and at the bottom of the domain for the B profile, while the thermal pressure of charges always remains below the magnetic pressure. The number densities (panel b) are different as compared to those of VALC due to the above mentioned reasons. However, given all the simplifications of our model, these differences are considered to be acceptable.

Fig. 1. Parameters of background model atmosphere as function of height. Panel a: gas pressure of neutrals (blue), charges (orange), and magnetic pressures corresponding to the S magnetic field profile (green) and B profile (red). Panel b: number density of neutrals (blue), and charges (orange). For comparison, corresponding number densities from VALC model atmosphere are plotted in green and red. Panel c: temperature. Panel d: horizontal magnetic field, Bx0, of the two background profiles, S and B. Here and below, zero of z axis is considered to be the same as in VALC atmosphere.

As is shown later, in order to understand the results of the simulations presented below, it is essential to know the characteristic frequencies associated with the background model atmosphere, such as the ion-cyclotron frequencies and the collision frequencies. The collision frequencies are calculated from Eq. (5), taking into account that ν_ni = ν_inn_e/n_n. These frequencies are given in Fig. 2 as a function of height. It is observed that for the S magnetic field profile, the ion-neutral collision frequency and the cyclotron frequency become equal in the lower part of the domain at z ≈ 0.8 Mm, while for the B magnetic field profile, the ion-cyclotron frequency is higher than the collision frequency throughout the domain. The highest wave frequency considered in our models corresponds to the period of 1 s (see below). This frequency is lower than any of the characteristic frequencies at all heights.

Fig. 2. Height dependence of characteristic frequencies calculated for background atmospheric model. The blue and orange lines are ion-cyclotron frequencies corresponding to the S and B magnetic field profiles, red and green lines are ion-neutral and neutral-ion collision frequencies, and violet is the highest wave frequency used is our work, corresponding to the wave period of 1 s.

2.3. Perturbation

The perturbation is generated by a driver at the base of the atmosphere. It consists of an oscillatory perturbation in the vertical fluid velocity and the corresponding perturbations in thermodynamic variables and magnetic field that are obtained from an analytical solution. Since the driver is located in the deep and dense layers, we assumed strong collisional coupling. The values we chose for the wave frequencies are much smaller than the collisional frequencies (see Fig. 2). Therefore, in order to generate the initial solution, we used single-fluid equations together with the generalized induction equation and kept only the leading effect due to ion-neutral interaction in the solar atmosphere, that is, the ambipolar diffusion,

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·\left(\rho \mathrm{u}\right)=0,\hfill \\ \hfill & \frac{\partial (\rho \mathit{u})}{\partial t}+\mathrm{\nabla }·(\rho \mathit{u}\mathit{u}+p)=\mathit{J}\times \mathit{B}+\rho \mathit{g},\hfill \\ \hfill & \frac{\partial }{\partial t}{\displaystyle (e+\frac{1}{2}\rho u^2)}+\mathrm{\nabla }·{\displaystyle (\mathit{u}(e+\frac{1}{2}\rho u^2)+p\mathrm{u})}=\mathit{J}·\mathit{E}+\rho \mathit{u}·\mathit{g},\hfill \\ \hfill & \frac{\partial \mathit{B}}{\partial t}-\mathrm{\nabla }\times (\mathit{u}\times \mathit{B})-\mathrm{\nabla }\times \left[{\eta }_A\frac{[(\mathit{J}\times \mathit{B})\times \mathit{B}]}{{|B|}^2}\right]=0,\hfill \end{array}$$(15)

where the ambipolar diffusivity coefficient is defined as

$$\begin{array}{c}\hfill {\eta }_A=\frac{{\xi }_{\mathrm{n}}^2{|B|}^2}{\alpha {\rho }_{\mathrm{n}}{\rho }_{\mathrm{c}}},\end{array}$$(16)

being ξ_n = ρ_n/ρ the neutral fraction, and the collisional parameter α as defined in Eq. (4).

These equations are linearized when taking into account that the magnetic field is purely horizontal and its variation with height is very small compared to the rest of background variables, so that the gradients in B_x0 can be neglected. It is also important to take into account that both the stratification and the perturbations are only functions of one vertical coordinate, z, and that another independent horizontal direction is x,

$$\begin{array}{cc}& \frac{\partial {\rho }_1}{\partial t}=-u_z\frac{\mathrm{d}{\rho }_0}{\mathrm{d}z}-{\rho }_0\frac{\partial u_z}{\partial z},\hfill \\ \hfill & {\rho }_0\frac{\partial u_z}{\partial t}=-{\rho }_{1}g-\frac{\partial p_1}{\partial z}-\frac{1}{{\mu }_0}\left(\frac{\partial B_{x1}}{\partial z}{B}_{x0}\right),\hfill \\ \hfill & \frac{\partial p_1}{\partial t}=c_{\mathrm{s}}^2\frac{\partial {\rho }_1}{\partial t}+c_{\mathrm{s}}^2{u}_z\frac{\mathrm{d}{\rho }_0}{\mathrm{d}z}-u_z\frac{\mathrm{d}{p}_0}{\mathrm{d}z},\hfill \\ \hfill & \frac{\partial B_{x1}}{\partial t}=-B_{x0}\frac{\partial u_z}{\partial z}+\frac{{\eta }_A}{{\mu }_0}\frac{{\partial }^2{B}_{x1}}{\partial z^2}+\frac{1}{{\mu }_0}\frac{\mathrm{d}{\eta }_A}{\mathrm{d}z}\frac{\partial B_{x1}}{\partial z},\hfill \end{array}$$(17)

where $c_{\mathrm{s}}=\sqrt{\gamma p_0/{\rho }_0}$ is the sound speed.

An analytical solution for these equations is still complicated and several simplifications need to be made. Since we needed to generate the analytical solution in only a few bottom grid points (the ghost points), we assumed that the temperature is locally uniform, that is, c_s is also locally uniform, and then the background pressure and density have the same vertical scale height H, as follows:

$$\begin{array}{c}\hfill \{p_0,{\rho }_0\}=\{p_0(z_0),{\rho }_0(z_0)\}\times exp(-z/H).\end{array}$$(18)

We further assumed that the ambipolar diffusion is negligible at the bottom layer. After these simplifications, Eq. (17) can be combined into a single wave equation, similar to Eq. (10) provided by Nye & Thomas (1976) for the particular case of k_x = k_y = 0,

$$\begin{array}{c}\hfill \frac{{\partial }^2{u}_z}{\partial t^2}=\frac{{\partial }^2{u}_z}{\partial z^2}(c_{\mathrm{s}}^2+{v_A}^2)-\frac{c_{\mathrm{s}}^2}{H}\frac{\partial u_z}{\partial z}\end{array}$$(19)

when calculating the solution in the few grid points at the bottom boundary, we assumed that the quantity $v_A\approx B_{x0}/\sqrt{{\mu }_0{\rho }_0}$ is uniform. Then, the dispersion relation is obtained as usual, assuming the solution of the form, u_z = Vexp(i(ωt − kz)),

$$\begin{array}{c}\hfill {\omega }^2=k^2(c_{\mathrm{s}}^2+{v_A}^2)-\frac{\mathit{ik}}{H}{c}_{\mathrm{s}}^2.\end{array}$$(20)

The polarization relations for the amplitudes of the rest of the variables can be derived from Eq. (17) by substituting,

$$\begin{array}{c}\hfill \{\frac{{\rho }_1}{{\rho }_0},\frac{p_1}{p_0},\frac{B_{x1}}{B_{x0}}\}=\{\stackrel{\sim }{R},\stackrel{\sim }{P},\stackrel{\sim }{B}\}exp(i(\omega t-kz)),\end{array}$$(21)

where $\{\stackrel{\sim }{R},\stackrel{\sim }{P},\stackrel{\sim }{B}\}$ are generally complex amplitudes. We obtained the following relations,

$$\begin{array}{cc}\hfill \stackrel{\sim }{R}& =\frac{V}{i\omega }{\displaystyle (ik+\frac{1}{H})},\hfill \\ \hfill \stackrel{\sim }{P}& =\frac{\gamma V}{i\omega }{\displaystyle (ik+\frac{1}{\gamma H})},\hfill \\ \hfill \stackrel{\sim }{B}& =\frac{\mathit{Vk}}{\omega }·\hfill \end{array}$$(22)

The wave period P = 2π/ω and the velocity amplitude V at the base of the atmosphere are free parameters of the simulations. Given ω, we obtained k from the dispersion relation, Eq. (20). We chose the velocity amplitude at the base of the atmosphere as a fraction of the background sound speed: V(z₀) = A × 10⁻³ ⋅ c_s(z₀), where the factor A was varied from 1 to 100. Also, the wave periods P = 1, 5, 7.5 and 20 s were used.

The initial condition for oscillations of neutrals and charges are set in the following way,

$$\begin{array}{cc}\hfill V_{\mathrm{n}}& =V;V_{\mathrm{c}}=V\hfill \\ \hfill {\stackrel{\sim }{R}}_{\mathrm{n}}& =\frac{{{\rho }_{\mathrm{n}}}_0}{{\rho }_0}\stackrel{\sim }{R};{\stackrel{\sim }{R}}_{\mathrm{c}}=\frac{{{\rho }_{\mathrm{c}}}_0}{{\rho }_0}\stackrel{\sim }{R};\hfill \\ \hfill {\stackrel{\sim }{P}}_{\mathrm{n}}& =\frac{{{\rho }_{\mathrm{n}}}_0}{{\rho }_0}\stackrel{\sim }{P};{\stackrel{\sim }{P}}_{\mathrm{c}}=\frac{{{\rho }_{\mathrm{c}}}_0}{{\rho }_0}\stackrel{\sim }{P}.\hfill \end{array}$$(23)

For our 1D numerical calculation, we covered the domain of L_z = 1.6 Mm with 32 000 grid points. The analytical wave solution was generated at each time step in the ghost points as the lower boundary condition. The perfectly matched layer (PML; see Berenger 1994; Felipe et al. 2010) was used as the upper boundary condition.

3. Nonlinear wave propagation in stratified plasma

The results of the numerical solution are given in Figs. 3 and 4. Throughout the simulations, we varied the following three parameters: the background magnetic field profile, the period, and the amplitude of the wave. The change in the background magnetic field and in the period of the wave can be studied in the linear approximation, while the change in the amplitude is related to the nonlinear effects. These figures show variations with the height of the charges, neutral velocities, their corresponding temperatures, and the magnetic field for a fixed moment of time. Figure 3 gives the results for two different initial wave amplitudes (factor A = 1 and 100) and two different periods (P = 20 and 7.5 s) for the background model with S magnetic field profile. Figure 4 shows oscillations for two magnetic field profiles (B and S) and two wave periods (P = 1 and 5 s) for the same wave amplitude factor A = 1.

Fig. 3. Height dependence of oscillations in velocity and temperature of charges and neutrals, and in magnetic field for fixed moment of time, obtained from numerical solution of two fluid equations for magnetic field profile S. Panels a and b: initial wave amplitude factor A = 1; panels c and d are for A = 100. The wave period is P = 20 s for panels a and c, and P = 7.5 s for panels b and d. Orange lines are for the parameters corresponding to charges, and blue lines are those for neutrals. We note that shock overshooting observed in panels c and d are numerical artifacts since we ran the simulation without any shock capturing algorithms in this case.

Fig. 4. Height dependence of oscillations in velocity and temperature of charges and neutrals, and in magnetic field for fixed moment in time, for initial wave amplitude factor A = 1. The format of the figure is the same as for Fig. 3. Panels a and b: S magnetic field profile; panels c and d are for the B profile. The wave period is P = 5 s for panels a and c, and P = 1 s for panels b and d.

The first impression from Fig. 3 is that the oscillation curves for neutrals and charges are nearly the same. The temperatures of both species are completely collisionally coupled via the thermal exchange process and, therefore, oscillations in both temperatures are the same. Nevertheless, one can observe some decoupling in the charges and neutral velocities at the upper part of the atmosphere, at heights above approximately 1.6 Mm. At these heights the collision frequency decreases because of the very low background densities. This strong drop of densities compensates for the increase in the background temperature and, consequently, in the value of the collisional parameter α. This decoupling is present, independently of the wave amplitude factor, and for both wave periods shown in the figure. It can be observed that velocity oscillations in neutrals slightly lag behind the oscillations in charges. The decoupling is visible not only at the shock wave fronts, but over the entire oscillation curve.

It is important to note that the vertical wavelength of the fast waves considered in this numerical experiment are dependent on the Alfvén speed. The increase of the Alfvén speed with height produces oscillations with significantly larger wavelength in the upper part of the atmosphere. This longer wavelength affects the formation of nonlinearities (panels c and d). Oscillations are strongly nonlinear in nature already at the bottom part of the domain for the case A = 100. Nevertheless, the shock fronts become smoother at bigger heights.

In general, velocity oscillations in Fig. 3 increase their amplitude with height. However, by comparing the panels a and c it can also be seen that the amplitude increase due to larger height is significantly lower, or even absent, in the case of strong perturbation. Perturbations in the temperature and magnetic field behave differently, and their amplitude starts to decrease roughly after the same height as where the charges-neutral velocity decoupling becomes visible.

When the wave period for the S magnetic field profile decreases further, a new effect becomes visible, namely oscillation damping (two upper panels of Fig. 4). As a result of this damping, the velocity amplitude decreases after about 1.4 Mm for P = 5 s wave, and above 1 Mm for P = 1 s wave. The damping is also apparently sensitive to the magnetic field profile. For the stronger field B (two bottom panels), the P = 1 s oscillation is damped from the very beginning, and the wave completely disappears after 0.8 Mm. For the P = 5 s, oscillations in the temperature and magnetic field also become visibly damped, while the velocity amplitude increase is still present. The difference between P = 1 and P = 5 is seen in significantly larger wavelength of oscillations in the case of the latter. This means that they are less affected by collisions. Furthermore, the effect on the damping is discussed in Sect. 4.

In summary, the following two main effects on waves are apparent from our simulations: decoupling in charges-neutral velocity and wave damping. These two effects are not independent from each other. The charges-neutral velocity decoupling can produce wave damping. However, the damping can also be produced by more mechanisms, which is demonstrated later in the paper. As for the case considered here, wave damping is a consequence of the charges-neutral decoupling and nonlinear effects, to a lesser extent.

3.1. Decoupling

The dependence of the degree of decoupling on the wave period and nonlinearities is difficult to evaluate just with the naked eye for Figs. 3 and 4. In order to quantify the effects, we adopted a measure of decoupling by computing the standard deviation of ∣u_cz − u_nz∣ in a stationary regime of the simulations (when the waves reached the upper boundary and several periods have passed through it), which is relative to the initial amplitude of the waves at the bottom level z₀. Such normalization gives us the possibility to better visualize the influence of the nonlinear effects. Figures 5–7 compare the amount of decoupling as a function of several parameters. In all of the cases discussed below, the amount of decoupling, normalized to u_cz(z₀), varies from a few percent of the initial wave amplitude at the lower layers to a significant fraction at the upper layers, which is due to the growth of the wave amplitude with height.

Fig. 5. Decoupling between charges and neutral velocity as function of height for magnetic field profile S. Left panel: wave period 20 s, amplitude factors A = 1 (blue), A = 10 (orange), A = 100 (green). Right panel: wave period 1 s, amplitude factors A = 0.5 (blue), A = 1 (orange), A = 2 (green).

Figure 5 reveals that the decoupling is a sensitive function of the wave amplitude. For the wave period of 20 s, when oscillation velocities are not visibly damped, the decoupling smoothly grows with height for all the wave amplitudes considered. For progressively more nonlinear cases (increasing A) the decoupling starts at lower heights and increases up to about 1.2 Mm. Above this height, nonlinear effects play an important role. Comparing the A = 1, A = 10, and A = 100 cases, the decoupling measured relative to the initial wave amplitude is smaller in the latter cases.

The action of nonlinear effects increases the relative difference between the velocity of charges and neutrals due to the formation of discontinuities associated with shocks. Nevertheless, a further increase of A reverses this tendency. The shock waves in the A = 100 case become strongly collisionally damped. Because of that, the amplitude of both u_cz and u_nz, and of their difference, decreases compared to the A = 10 case. The values of the decoupling measured with respect to the initial wave amplitude become lower at heights where strong shock damping takes place. However, the absolute value of the decoupling, without normalization, are still larger in the strongly nonlinear case A = 100.

The right panel of Fig. 5 shows the case of P = 1 s where strong wave damping was observed. In this case, the decoupling has a broad maximum between 0.8 and 1.2 Mm, and is zero approximately above 1.4 Mm since no waves exist at higher heights. Similar to the previous case, the amount of decoupling depends on the wave amplitude. We observed that the peak becomes slightly lower and it is progressively shifted to lower heights with an increasing A due to nonlinear effects.

The amount of decoupling also depends on the wave period. As is seen in Fig. 6, the values of the decoupling smoothly decrease with the period for P = 5−20 s, showing the same height dependence. The decoupling is maximum for the strongly damped case, P  =  1 s at intermediate heights; however, it disappears completely above a certain height for the reasons explained above. The increase of the magnetic field causes lower decoupling, as shown in Fig. 7, and shifts the peak to lower heights.

Fig. 6. Decoupling between charges and neutral velocity as function of height for magnetic field profile S and wave amplitude factor A = 1. Wave periods P = 20 s (violet), P = 7.5 s (red), P = 5 s (green), P = 2.5 s (orange), and P = 1 s (blue).

Fig. 7. Decoupling between charges and neutral velocity as function of height for amplitude factor A = 1 and wave period P = 1 s. The dependence on the magnetic field profiles is shown with different curves: S (blue) and B (orange).

To better understand the role of the nonlinear effects and collisions for the decoupling, we ran the simulation with P  =  1 s again. This time we only evolved linear equations with either the original collisional parameter α (Eq. (4)), corresponding to the atmospheric model, or artificially increased the collisions by a factor of 10⁴. The results of this experiment are given in Fig. 8, together with the original nonlinear case. It can be seen that the decoupling is larger and its maximum is shifted to higher heights in the linear case (orange line), showing the same tendency as in Fig. 5 (right panel). In the case when the collisional parameter α is artificially increased, the decoupling completely disappears (green line).

Fig. 8. Decoupling between charges and neutral velocity as function of height for S magnetic field profile and amplitude factor A = 1. Blue line: nonlinear simulation; orange line: linear simulation; green line: linear simulation with artificially increased collisional parameter α.

Therefore, the overall conclusion from the analysis above is that the decoupling increases with decreasing the wave period due to the relative increase of the importance of the collisional time scales compared to the wave scales. Similarly, the absolute value of the decoupling increases with the wave amplitude since nonlinearities bring shorter scales to the problem due to the formation of shock waves. However, the value of the decoupling relative to the initial wave amplitude is lower in the strongly nonlinear cases due to the collisional damping of the shock waves and the subsequent decrease of the oscillation amplitude of neutrals and charges.

3.2. Evolution of background variables

The difference in the velocity of charges and neutrals leads to frictional dissipation of the kinetic energy of both fluids and, therefore, can produce net variations in the atmospheric parameters, such as the temperature. To check the amount of these net variations, we computed the time averages of the perturbed variables as a function of height. The results are given in Fig. 9 for the case of P = 1 s, where the maximum damping was observed for the S magnetic field profile.

Figure 9 shows that the temperature has, on average, increased at heights between 0.8 and 1.2 Mm, which coincides with the height range where the decoupling was maximum for this case (Fig. 5). The amount of the temperature increase is about 0.5 K in this case; however, it must be compared to the amplitude of the temperature oscillations (about 10 K). Therefore, the temperature increase is not negligible. The average temperature increase is accompanied by the decrease of the magnetic field strength and variations of density and net velocity of the charges and neutrals. As the temperature increases with time, the pressure increases as well producing gas expansion and a corresponding decrease in the densities. The plasma tied to the magnetic field expands and compresses the field lines. Therefore, the increase in temperature is accompanied by the decrease in the magnetic field. The existence of the non-zero net flow in neutrals and charges (bottom panel) confirms this picture.

Fig. 9. Time average of perturbed variables as function of height for simulation with S magnetic field profile, wave period P = 1 s and amplitude factor A = 2. From top to bottom: temperatures of neutrals (orange) and charges (blue), magnetic field, density of charges, density of neutrals, velocity of neutrals (orange) and charges (blue).

A better understanding of the above result is gained through considering the equations of evolution of the internal energy, which are obtained after removing the kinetic energy part from the corresponding equations of the two-fluid system (Eq. (1)),

$$\begin{array}{cc}\hfill \frac{\partial e_{\mathrm{n}}}{\partial t}& +\mathrm{\nabla }·\left({\mathit{u}}_{\mathrm{n}}{e}_{\mathrm{n}}\right)+p_{\mathrm{n}}\mathrm{\nabla }·{\mathit{u}}_{\mathrm{n}}=Q_{\mathrm{n}},\hfill \\ \hfill \frac{\partial e_{\mathrm{c}}}{\partial t}& +\mathrm{\nabla }·\left({\mathit{u}}_{\mathrm{c}}{e}_{\mathrm{c}}\right)+p_{\mathrm{c}}\mathrm{\nabla }·{\mathit{u}}_{\mathrm{c}}=\mathit{J}·{\mathit{E}}^{\ast }+Q_{\mathrm{c}},\hfill \end{array}$$(24)

where E^* = [E + u_c × B]. The expressions for collisional terms, Q_n and Q_c have the following form,

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

These collisional terms can be compared with M_n in the total energy equation. The Q_n and Q_c terms do not strictly compensate each other as is the case of M_n. There are contributions proportional to the square of the velocity difference, (u_c − u_n)², that are positive and are added to both charges and neutral energy equations, called frictional heating (FH). The second term, that is, the thermal exchange (TE), has the same absolute value, but different sign for the neutrals and charges.

In order to translate the contribution of Q-terms into the temperature increase, one can rewrite the internal energy equation in terms of temperature by substituting e_n = n_nk_bT/(γ − 1) and e_c = 2n_ek_bT/(γ − 1) = n_ck_bT/(γ − 1). By operating the energy equation using the continuity equations, the temperature evolution equations are obtained,

$$\begin{array}{cc}\hfill \frac{\partial T_{\mathrm{n}}}{\partial t}+\mathrm{\nabla }·\left({\mathit{u}}_{\mathrm{n}}{T}_{\mathrm{n}}\right)& +(\gamma -1)T_{\mathrm{n}}\mathrm{\nabla }·{\mathit{u}}_{\mathrm{n}}=Q_{\mathrm{n}}\frac{\gamma -1}{k_b{n}_{\mathrm{n}}},\hfill \\ \hfill \frac{\partial T_{\mathrm{c}}}{\partial t}+\mathrm{\nabla }·\left({\mathit{u}}_{\mathrm{c}}{T}_{\mathrm{c}}\right)& +(\gamma -1)T_{\mathrm{c}}\mathrm{\nabla }·{\mathit{u}}_{\mathrm{c}}=(\mathit{J}·{\mathit{E}}^{\ast }+Q_{\mathrm{c}})\frac{\gamma -1}{k_b{n}_{\mathrm{c}}}·\hfill \end{array}$$(26)

The contributions to the terms Q_n, c(γ − 1)/k_bn_n, c are computed from the same simulation as above. These terms appear to be constant in time. Figure 10 shows the height dependence of the frictional heating and thermal exchange terms from Eq. (26), for neutrals (upper panel) and charges (middle panel), together with their joint contribution (bottom panel). As expected, the thermal exchange is a positive quantity for the neutrals and a negative quantity for the charges. The negative sign means that the energy would be transferred from charges to neutrals. Due to significantly lower number density of charges, the absolute values of both FH and TE terms in the temperature equation for charges are three orders of magnitude larger than for neutrals. Therefore, the collisions lead to efficient frictional heating of charges, and then this energy is transferred to neutrals via the thermal exchange. With the charges and neutrals thermally coupled, the evolution of the temperature of the fluid as a whole is dominated by the evolution of the neutral temperature.

Fig. 10. Dissipative terms in temperature equations (Eq. (26)) as function of height for simulation with S magnetic field profile, wave period P = 1 s, and amplitude factor A = 2. Top panel: thermal exchange term (TE, orange line) and frictional heating term (FH, blue line) for neutrals. Middle panel: same for charges. Bottom panel: total energy exchange terms for charges (blue), neutrals (orange). We note that apparent oscillations in the bottom panel are numerical artifacts due to averaging the signal with relatively low temporal resolution.

The constant action of the frictional heating terms leads to an uninterrupted energy supply and an increase of temperature in time. Figure 11 shows the time evolution of temperature (left panel) and the magnetic field (right panel) at a fixed point z = 1 Mm for the case P = 1 s, S magnetic field profile, and varying the amplitude of the perturbation, A. The curves are obtained through calculating the running average over three consecutive wave periods. It can be observed that the increase in temperature is a time cumulative effect, proportional to the square velocity amplitude difference, and that the temperature grows linearly in time. The growth is linear while the modification of the background is of the same order or smaller than the perturbation amplitude. This linear growth is consistent with the value of Q_n(γ − 1)/k_bn_n ≈ 8 × 10⁻³ K s⁻¹, see the bottom panel of Fig. 10 at the location z = 1 Mm. The corresponding drop of the magnetic field strength is also linear with time. The gas is heated and it expands, expelling the field lines, which is consistent with the picture from Fig. 9.

Fig. 11. Running average of perturbed temperature (left) and magnetic field (right) as function of time, for fixed height z = 1 Mm for simulation with S magnetic field profile, wave period P = 1 s. The averaging is taken over three consecutive wave periods. The amplitude factor is A = 0.5 (blue), A = 1 (green), and A = 2 (red). Temperature of neutrals and charges closely match together.

4. Reasons for wave damping

One of the interesting results obtained above is the wave damping, that is, the decrease of the wave amplitude (or even a complete disappearance of the wave) after a certain height. This effect is especially pronounced for the high-frequency perturbations.

There can be several reasons for this damping. One of them can be ion-neutral collisions, which are well investigated in linear theory (Ballester et al. 2018a). Another reason can be wave reflection due to the gradient of the Alfvén speed with height. It is well known that fast magneto-acoustic waves refract and reflect in the solar atmosphere due to this reason (Cally 2006; Khomenko & Collados 2006). This is a phenomenon that is actually thought to be responsible for the formation of high-frequency acoustic halos in observations (Khomenko & Collados 2009; Rijs et al. 2016).

In order to check if the waves are reflected in our experiment, we computed the acoustic wave energy fluxes, as p_c1u_zc for charges, and as p_n1u_zn for neutrals. It turns out that both fluxes are positive at all heights, meaning vertical wave propagation and no apparent reflection. The energy fluxes vanish after a certain height above 1.2 Mm, which means that the waves became non-propagating. The reasons of this behavior are better highlighted with the help of the analytical solution, developed below.

In order to distinguish between the damping produced by the decoupling of the neutral and charged fluids in the linear regime, and the damping caused by the decoupling at the shock fronts or the steepening in the wave profile in the nonlinear regime, we compared the numerical solution of the linear and nonlinear two-fluid equations with the analytical solution. The analytical solution is obtained by solving the linearized two fluid equations with a nonuniform background for each species. We introduced the coupling term in the momentum equations and assumed the collisional parameter α to be constant in time, which is equal to the value corresponding to the equilibrium atmosphere and consistent with the linear approximation. The linearized equation takes the following form,

$$\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}}{\partial z}-u_{\mathrm{cz}}\frac{\mathrm{d}{p}_{\mathrm{c}0}}{\mathrm{d}z},\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}}{\partial z}-u_{\mathrm{nz}}\frac{\mathrm{d}{p}_{\mathrm{n}0}}{\mathrm{d}z},\hfill \\ \hfill & \frac{\partial B_{x1}}{\partial t}=-B_{x0}\frac{\partial u_{\mathrm{cz}}}{\partial z}-u_{\mathrm{cz}}\frac{\mathrm{d}{B}_{\mathrm{x}0}}{\mathrm{d}z}·\hfill \end{array}$$(27)

These equations are combined into two coupled partial differential equations for the vertical velocity of each species and coupled by the collisional terms as follows:

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

where

$$\begin{array}{cc}& a_{\mathrm{c}}(z)=c_{\mathrm{c}0}^2+{v_A}_0^2,a_{\mathrm{n}}(z)=c_{\mathrm{n}0}^2,\hfill \\ \hfill & b_{\mathrm{c},\mathrm{n}}(z)=\frac{1}{{{\rho }_{\mathrm{c},\mathrm{n}}}_0}\frac{\mathrm{d}\left({{\rho }_{\mathrm{c},\mathrm{n}}}_0{a}_{\mathrm{c},\mathrm{n}}\right)}{\mathrm{d}z},\hfill \\ \hfill & {c_{\mathrm{c},\mathrm{n}}}_0^2=\gamma \frac{{p_{\mathrm{c},\mathrm{n}}}_0}{{{\rho }_{\mathrm{c},\mathrm{n}}}_0},{v_A}_0^2=\frac{{B_x}_0^2}{{\mu }_0{{\rho }_{\mathrm{c}}}_0}·\hfill \end{array}$$(29)

Assuming the functional dependence on space and time of the vertical velocity of the form $\{u_{\mathrm{cz}}(z,t),u_{\mathrm{nz}}(z,t)\}=\{{\stackrel{\sim }{u}}_{\mathrm{cz}}(z),{\stackrel{\sim }{u}}_{\mathrm{nz}}(z)\}\times \text{exp}\left(i\omega t\right)$, where the frequency ω is fixed, the above coupled system can be rewritten as,

$$\begin{array}{cc}\hfill -{\omega }^2{\stackrel{\sim }{u}}_{\mathrm{cz}}& =a_{\mathrm{c}}(z)\frac{{\mathrm{d}}^2{\stackrel{\sim }{u}}_{\mathrm{cz}}}{\mathrm{d}{z}^2}+b_{\mathrm{c}}(z)\frac{\mathrm{d}{\stackrel{\sim }{u}}_{\mathrm{cz}}}{\mathrm{d}z}+i\omega \alpha {\rho }_{\mathrm{n}0}({\stackrel{\sim }{u}}_{\mathrm{nz}}-{\stackrel{\sim }{u}}_{\mathrm{cz}}),\hfill \\ \hfill -{\omega }^2{\stackrel{\sim }{u}}_{\mathrm{nz}}& =a_{\mathrm{n}}(z)\frac{{\mathrm{d}}^2{\stackrel{\sim }{u}}_{\mathrm{nz}}}{\mathrm{d}{z}^2}+b_{\mathrm{n}}(z)\frac{\mathrm{d}{\stackrel{\sim }{u}}_{\mathrm{nz}}}{\mathrm{d}z}+i\omega \alpha {\rho }_{\mathrm{c}0}({\stackrel{\sim }{u}}_{\mathrm{cz}}-{\stackrel{\sim }{u}}_{\mathrm{nz}}).\hfill \end{array}$$(30)

In combining these two equations, we obtained the folllowing:

$$\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{cz}}}{\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 }_{n0}))\hfill \\ \hfill & +{\stackrel{\sim }{u}}_{\mathrm{cz}}{\omega }^3(\omega -i\alpha ({\rho }_{\mathrm{c}0}+{\rho }_{\mathrm{n}0}))=0.\hfill \end{array}$$(31)

The coefficients of the spatial derivatives that appear in Eq. (31) are not uniform and there is no apparent exact analytical solution to Eq. (31). The waves we consider in this study are short-period waves with wavelengths that are shorter than equilibrium gradient scales. Thus, we used the WKB approximation to find solutions for Eq. (31), as follows:

$$\begin{array}{c}\hfill \{{\stackrel{\sim }{u}}_{\mathrm{nz}},{\stackrel{\sim }{u}}_{\mathrm{cz}}\}=\{V_{\mathrm{n}}(z),V_{\mathrm{c}}(z)\}·\text{exp}[i\varphi (z)],\end{array}$$(32)

where V_n and V_c are space dependent complex amplitudes, and

$$\begin{array}{c}\hfill \varphi (z)=-{{\displaystyle \int }}_0^{z}k(z^{\prime })\mathrm{d}{z}^{\prime }.\end{array}$$(33)

Similar to the WKB approximation, we assumed that the velocity amplitude and the wavenumber gradients are small and of the same order,

$$\begin{array}{cc}& \frac{\mathrm{d}\varphi }{\mathrm{d}z}=-k=-k_0-ϵk_1(z),\hfill \\ \hfill & V_{\mathrm{c},\mathrm{n}}(z)=V_{\mathrm{c},\mathrm{n}0}+ϵV_{\mathrm{c},\mathrm{n}1}(z),\hfill \end{array}$$(34)

with ϵ being a small parameter and the second derivatives of k₁(z) and V_c, n1 being zero. It is then straightforward when making the following calculations:

$$\begin{array}{cc}& \frac{\mathrm{d}{\stackrel{\sim }{u}}_{\mathrm{cz}}}{\mathrm{d}z}=(ϵ\frac{\mathrm{d}{V}_{c1}}{\mathrm{d}z}-ik{V}_{c1})·\text{exp}(i\varphi ),\hfill \\ \hfill & \frac{{\mathrm{d}}^2{\stackrel{\sim }{u}}_{\mathrm{cz}}}{\mathrm{d}{z}^2}=(-2iϵk\frac{\mathrm{d}{V}_{c1}}{\mathrm{d}z}+iϵV_{c1}\frac{\mathrm{d}{k}_1}{\mathrm{d}z}-k^2{V}_{\mathrm{c}1})·\text{exp}(i\varphi ),\hfill \\ \hfill & \frac{{\mathrm{d}}^3{\stackrel{\sim }{u}}_{\mathrm{cz}}}{\mathrm{d}{z}^3}=(-3ϵk\frac{\mathrm{d}{V}_{c1}}{\mathrm{d}z}-3ϵV_{c1}\frac{\mathrm{d}{k}_1}{\mathrm{d}z}+i{k}^2{V}_{\mathrm{c}1})·k·\text{exp}(i\varphi ),\hfill \\ \hfill & \frac{{\mathrm{d}}^4{\stackrel{\sim }{u}}_{\mathrm{cz}}}{\mathrm{d}{z}^4}=(4iϵk\frac{\mathrm{d}{V}_{c1}}{\mathrm{d}z}+6iϵV_{c1}\frac{\mathrm{d}{k}_1}{\mathrm{d}z}+k^2{V}_{\mathrm{c}1})·k^2·\text{exp}(i\varphi ).\hfill \end{array}$$(35)

By substituting these expressions into Eq. (31) and in only keeping 0th-order terms, the dispersion relation for k is obtained,

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

where the remaining coefficients are defined as,

$$\begin{array}{c}\hfill b=({\rho }_{\mathrm{n}0}{b}_{\mathrm{n}}+{\rho }_{\mathrm{c}0}{b}_{\mathrm{c}})/{\rho }_0,a=({\rho }_{\mathrm{n}0}{a}_{\mathrm{n}}+{\rho }_{\mathrm{c}0}{a}_{\mathrm{c}})/{\rho }_0\end{array}$$(37)

with ρ₀ = ρ_c0 + ρ_n0.

Next, the relations (35) are replaced into Eq. (30). Separating the 0th-order terms from the 1st-order terms, the expressions for the velocity amplitudes are obtained as follows:

$$\begin{array}{c}\hfill V_{\mathrm{c},\mathrm{n}}(z)=V_{\mathrm{c},\mathrm{n}}{|}_{z=0}·\text{exp}\left({{\displaystyle \int }}_0^z\frac{\mathrm{d}k}{\mathrm{d}z}\frac{i{a}_{\mathrm{c},\mathrm{n}}}{b_{\mathrm{c},\mathrm{n}}-2ik{a}_{\mathrm{c},\mathrm{n}}}\mathrm{d}{z}^{\prime }\right).\end{array}$$(38)

In the limit of $\frac{\alpha {\rho }_0}{\omega }\gg 1$, Eq. (36) reduces to the single-fluid dispersion relation,

$$\begin{array}{c}\hfill -{\omega }^2+k^{2}a+ikb=0.\end{array}$$(39)

In the opposite limit, when $\frac{\alpha {\rho }_0}{\omega }\ll 1,$ we recovered the dispersion relations of either neutrals or charges. Otherwise, if $\frac{\alpha {\rho }_0}{\omega }\approx 1$, both of the terms in the equation are equally important.

It is interesting to analyze the effect of collisional damping separately from the effect of stratification. The effect of stratification is independent of the wave frequency for high-frequency waves considered here. By neglecting the stratification, the dispersion relation, Eq. (36) becomes the dispersion relation in a uniform atmosphere,

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

We considered the solutions of this equation as a function of the wave period. Using the non-dimensional variables,

$$\begin{array}{c}\hfill F=\alpha {\rho }_0/\omega ;E=k\sqrt{a}/\omega ,\end{array}$$(41)

the dispersion relation, Eq. (40), becomes

$$\begin{array}{c}\hfill (1-E^2{a}_{\mathrm{c}}/a)(1-E^2{a}_{\mathrm{n}}/a)-iF(1-E^2)=0.\end{array}$$(42)

Figure 12 shows the real and imaginary part of the four solutions, E, of the fourth order dispersion relation, Eq. (42), as a function of the non-dimensional variable F, which is the ratio between the collisional and the wave frequencies, scaled with ρ₀. This dispersion relation was solved using the values for a_n, a_c, and a corresponding to a point located in the middle of the atmosphere for S and B magnetic field profiles. The values on the horizontal axis cover a large range, from F = 10⁻⁴, for the uncoupled case, and up to F = 10⁵, for the strongly coupled case. We marked the values of F corresponding to the cases of ω = ν_in (ion-neutral collision frequency), ω = ν_ni (neutral-ion collision frequency), and ω = 2π/1 (frequency corresponding to the period of a 1 s wave) with F_c, F_n, and F₁, respectively. The values of the frequencies were calculated for the height corresponding to the middle of the domain.

We can observe that for small F, corresponding to high frequencies ω (uncoupled case), the waves propagate with the characteristic speed of either charges ($\sqrt{a_{\mathrm{c}}}$, downward solution #1, and upward solution #4) or with the characteristic speed of neutrals ($\sqrt{a_{\mathrm{n}}}$, downward solution #2, and upward solution #3). The imaginary part of k in the uncoupled case is zero, meaning that there is no damping. This result is similar to Popescu Braileanu et al. (2019).

For the large values of F, corresponding to low frequencies ω (strongly coupled case), there are two solutions, which propagate with the characteristic speed in the fluid as a whole, $\sqrt{a}$, both upward ( #2) and downward (#3). For the other two solutions, #1 and #4, the real and the imaginary parts of k are very large. These are waves with very large phase speed, but also extremely large damping, and we consider them as unphysical.

Solutions #2 and #3 (yellow and green lines) have maximum damping, k_I, for waves in the frequency range between ν_ni and ν_in. This is true for both magnetic field profiles B and S. For the larger magnetic field B, the maximum damping shifts to higher frequencies. We checked that for even higher magnetic fields the dispersion curves remain practically equal to the case B. For the 1 s period wave, the value of damping is decreased from its maximum, but it is still non-zero, which explains the damping we observe in our simulations.

The full dispersion relation, which includes stratification, (36) is also fourth order in k and has four different roots. Examples of these four different solutions as a function of the height in the atmosphere are illustrated in Fig. 13 and were calculated for the S magnetic field profile and the wave period of 5 s. Accordingly to Fig. 12, from the point of view of the collisions, this is the coupled case. It can be seen that the solutions #1 and #4 both have very large imaginary and real parts. These modes therefore have very low propagation speeds and very high damping, which is similar to the case without stratification. We did not consider these solutions further. The other two solutions marked #2 and #3 represent waves that propagate either up (solution #3) or down (solution #2) and their damping is moderate. It can be seen that k_I of both solutions matches for the first 0.5 Mm from the bottom, and has positive values, corresponding to the wave amplification. At these heights, the charges-neutral decoupling is not strong enough to damp the waves and the growth of the wave amplitude with height is caused simply by the gravitational stratification of the atmosphere. Higher up, k_I turns negative for the solution #3, meaning wave damping. The same amount of damping is also present in the solution #2, however k_I is positive because the wave propagates downward (negative k_R). For the comparison below, we use the solution #3.

Fig. 12. Solutions of non-dimensional dispersion relation without stratification, Eq. (42). Upper and lower panels: real and imaginary part of k, scaled with $\sqrt{a}/\omega$ as a function of the ratio between the collisional and wave frequencies, scaled with ρ0. Solutions labeled as “S” (dashed lines) are for S magnetic field; solutions labeled “B” (solid lines) are for B magnetic field. The vertical dotted lines mark values of ω = νin (Fc), ω = νni (Fn), and ω = 2π/1 (F1) for the 1 s period wave.

Fig. 13. Solutions of dispersion relation (36) for two-fluid linearized equations. Upper and lower panels: real and imaginary part of k as a function of height, calculated for the magnetic field profile S and 5 s period.

Figure 14 compares the numerical solutions of the linearized equations with the analytical solution for different wave periods and the magnetic field profile S. One can observe that both solutions match rather well and small differences are only observed in the upper layers. The analytical solution in the WKB approximation is able to describe rather precisely the damping of waves due to collisions in the linear regime. Therefore, we can conclude that the wave damping observed in these simulations is determined mainly by the linear decoupling.

Fig. 14. Comparison between numerical solution in linear regime (black lines) and analytical solution of two-fluid equations (red dotted line). Individual panels show snapshots of the velocity of charges as a function of height at fixed time moments in the stationary regime of the simulations. Panels from left to right, from top to bottom are for wave periods of 1, 5, 7.5, and 20 s, and S magnetic field profile. It is important to note that the analytical solution slightly overshoots the numerical solution, the effect being more pronounced toward the upper part of the domain. This is due to the fact that the analytical solution is an approximate solution.

Figure 15 further illustrates the effects of linear decoupling on the wave amplitudes. It shows the imaginary part of the wave vector, k_I, for the solution #3 as a function of the wave period (left) and of the magnetic field strength (right). The value and sign of k_I is directly related to the modification of the wave amplitude with height. It can be seen that in all cases, k_I changes sign at some height in the atmosphere. At the bottom part of the atmosphere, k_I is positive, meaning that the growth of the wave amplitude with height is due to gravitational stratification, as explained above. Once the collision frequency decreases with height, the effects of the linear wave damping become progressively more important and k_I changes sign. The damping starts at lower heights for shorter period waves. This happens because the wave frequency approaches the neutral-ion collisional frequency (ν_ni) and the collisions alter the wave propagation, consistently with Fig. 12, which is where we did not take the stratification into account. The damping also starts at lower height for stronger magnetic fields, which is consistent with the results of the simulations shown in the previous sections. However, the absolute value of damping is higher for lower B (so the damping length is shorter). Decreasing the value of the magnetic field for a fixed wave frequency results in a decrease of the wavelength, which facilitates the effects of collisions on the wave damping.

Fig. 15. Left panel: imaginary part of k as function of wave period (horizontal axis) and height (vertical axis) obtained after solving dispersion relation for two-fluid equations (Eq. (36)) for S magnetic field profile. Solution #3 from Fig. 13 is used. The period varies between 1 s and 20 s (the minimum and maximum values used in the simulations). The solid, black line represents the contour of kI = 0. Right panel: imaginary part of k as function of background magnetic field at base of atmosphere (horizontal axis) and height (vertical axis). The wave period is fixed to 5 s.

Considering the propagation for different magnetic field strengths, we observed the following two effects: one due to the stratification, which increases the amplitude of the wave, and the other due to the collisions which damps the wave. The amplification effects are stronger in the lower part of the atmosphere because of the lower temperatures and, consequently, the pressure scale height. The damping effects are stronger in the upper part of the atmosphere. The interplay between two effects can be observed in Fig. 13 by analyzing the imaginary part of solutions #2 and #3. In the lower part of the atmosphere the values are positive (meaning an increase of the amplitude) because of the stratification, and are identical for #3 and #2 because the effects of collisional damping are negligible. When the collisional damping become important in the upper part of the atmosphere, the two curves split and become symmetric with respect to the k_I = 0 line, meaning that the stratification does not play a major role in amplifying the wave amplitude anymore. The absolute values of both, amplification and damping are smaller for larger magnetic fields. The height when the damping effects overcome those of the stratification is lower for smaller magnetic fields.

In approaching the limit of the absence of the magnetic field, B ≈ 0 since the temperature of charges and neutrals are the same. The propagation speeds of both species also become very close to each other. In this case, there is very little damping for the periods considered here, as we can see in panel b of Fig. 15.

The nonlinear effects on the wave damping are shown in Fig. 16 . This figure compares the numerical solutions done with the same initial amplitude (A = 10), but in the linear and nonlinear regimes. The nonlinear effects significantly increase the wave damping for all considered wave frequencies. At lower heights, when the nonlinear steepening of the wave profile is not yet pronounced, the linear and nonlinear solutions match, both showing an increase of the wave amplitude with height due to the gravitational stratification. However, at higher heights, the amplitudes of the linear and nonlinear solutions become very different. The waves are significantly more damped in the nonlinear regime, and the effect strongly depends on the wave period. The effective shortening of the scale at the wave fronts facilitates the collisional damping. The nonlinear fronts form at lower heights for shorter period waves and, therefore, nonlinear effects for such waves are significantly more pronounced. One can compare the effective decrease of the amplitude between the nonlinear and linear cases in Fig. 16, showing how the progressively larger periods are less affected by the nonlinear effects. We conclude that nonlinear steepening of the wave fronts dramatically increases the collisional damping of waves in the chromosphere.

Fig. 16. Comparison between numerical solution in linear regime (black lines) and nonlinear regime (red dotted line) with A = 10. Individual panels show snapshots of the velocity of charges as a function of height at fixed time moments in the stationary regime of the simulations. Panels from left to right, from top to bottom are for wave periods of 1, 5, 7.5, and 20 s, and S magnetic field profile.

5. Discussion and conclusions

In this work we performed simulations of fast magneto-acoustic waves in the solar chromosphere using a two-fluid approach. We used stratified temperature distribution and number densities for charges and neutrals consistent with VALC solar model. Additionally, we considered two magnetic field profiles with an order of magnitude of different field strengths. Our main findings can be summarized as follows,

• Waves in neutrals and charges, initially coupled in the photosphere, become decoupled at some height in the chromosphere. This means that velocities of co-located fluid elements for neutrals and charges are different. In our simulations, the decoupling happens above 0.7–1 Mm height. The difference in the charges-neutral velocity can reach up to 30–50% of the initial wave amplitude at the bottom of the domain (about 0.5 Mm in our model).

• The temperatures of neutrals and charges are effectively coupled by the thermal exchange and oscillate with the same amplitude and phase.

• The decoupling in the wave velocity is a function of the wave period, its amplitude, and the background magnetic field strength. In general, waves with smaller periods show greater decoupling at the same height. Decoupling happens at lower heights for stronger magnetic fields.

• Charge-neutral collisions cause significant wave damping. We observe significant damping of the shorter period waves, 1 and 5 s in our simulations. In such cases, the absolute value of the decoupling is smaller compared to the less damped case because the velocity amplitude strongly decreases after some height in the chromosphere and, in some cases, the perturbation associated with the waves completely disappears. Waves are damped at lower heights when the magnetic field is larger, but the damping length is shorter in the case of the smaller magnetic field.

• The damping obtained in the simulations in the linear regime can be rather precisely described by the analytic theory solving linear two-fluid equations using the WKB approach. This level of agreement suggests that the damping observed in such simulations is determined by the linear effects due to decoupling.

• Nonlinear wave propagation effects and steepening of wave fronts largely enhance the effects of the linear collisional damping. Nonlinear effects for the damping are more pronounced for the shorter period waves because the shock formation happens at lower heights.

• Collisional damping produces an effective frictional heating and local background temperature increase. The rate of the temperature increase is constant in time and is proportional to the square of the velocity amplitude.

The charges-neutral velocity decoupling and wave damping are two related multi-fluid effects. The decoupling appears when the collisional time scales become similar or longer than the hydrodynamic time scales. The decoupling determines damping of the waves. The kinetic energy lost in the work done by the collisional terms is converted into internal energy through the frictional heating term, which is a positive quantity added in the same amount to the energy equation of both species. The damping of several types of waves by collisions is also suggested in other studies (e.g., Zaqarashvili et al. 2011, 2013; Soler et al. 2013a,b).

In our case, the hydrodynamic time scales for the charges are determined by the ion cyclotron frequency and by the wave frequency. The wave frequencies considered here are much smaller than the ion cyclotron frequency at all heights. The collision time scales are determined by the ion-neutral and neutral-ion collision frequencies for the charges and neutrals, respectively. The relative importance of collisional terms, which couple the evolution of charges and neutrals, varies strongly with height due to the gravitational density stratification, as shown in Fig. 2. For the small background magnetic field, S profile, the ion-neutral collision frequency becomes smaller than the cyclotron frequency at z  ≈  0.8 Mm in the atmosphere, and we expected the ions to decouple from the neutrals above this point, leading to wave damping. With the same argument, but for the larger magnetic field (B profile), we expected the waves to be damped at even lower heights than for the S profile because the ion cyclotron frequency increases and equals the collision frequency at the bottom of the atmosphere. Our numerical simulations show that this is indeed the case. However, there is an additional effect caused by the gravitational stratification of the atmosphere, which can produce strong variation of the perturbation wavelength with height, depending on the magnetic field strength. The perturbation wavelength is larger for larger fields. Consequently, since the perturbation gradient spatial scale is larger, the absolute value of the damping in the strong-field case is smaller (or, in other words, the damping length is larger). Both numerical simulations and the analytical WKB solution show similar values of the damping.

In the case of the stratified atmosphere, the effect of the damping competes with the effect of the wave amplitude increase due to the gravitational stratification. The exponential decrease of the background density with height results in the increase of the wave amplitude at lower heights. This effect is well reproduced in the analytical WKB solution.

The nonlinear effects result in the steepening of the wave fronts and in the formation of shocks. The effective spatial scale at shock fronts decreases and, depending on the wave period and amplitude, it can reduce to the collision mean free path between neutrals and charges, which is larger in the upper part of the atmosphere, and attains the value of 2 km at the height of 1.5 Mm. But even if the shock width is several orders of magnitude larger than the mean free path, as in the case of Hillier et al. (2016), its size, that is much smaller than the wavelength, greatly enhances the damping of the waves. This effect is observed in the simulations when comparing the numerical solutions in the linear and nonlinear regimes in Fig. 16.

If the decrease in the amplitude of the shock due to collisional damping is large enough, the discontinuity at the shock front smoothes. This smoothing can be observed in Fig. 16 (upper right and bottom left panels). The wave fronts have a marked saw-tooth profile at intermediate heights, but they are smooth at the upper part of the atmosphere. These smoothed shock profiles resemble a C-type shock. Shock waves in a weakly collisional interstellar medium were extensively studied by Mullan (1971), Draine (1980), and Draine et al. (1983). It was found that when the Alfvén speed is larger than the shock speed, these shocks are preceded by a “magnetic precursor”, which heats and compresses the medium ahead of the front where the neutral gas undergoes a discontinuous change of state. If the magnetic field is strong enough, the shock will be a C-type shock (when all the variables are continuous across the shock). Otherwise, the neutral gas variables are discontinuous (J-type shock). A substantial fraction of the energy would be dissipated in this magnetic precursor because of the ion-neutral streaming. The solar atmosphere differs from the interstellar medium due to significantly higher density and temperature; therefore, the decoupling is less pronounced. Nevertheless, Hillier et al. (2016) find the existence of two both C-type and J-type shocks in their study of the reconnection-driven, slow-mode shocks in the solar corona. Hillier et al. (2016) observe a continuous transition from C-type shocks for subsonic velocity upstream of the shock to J-type shocks, for all the variables, in the case of supersonic upstream velocity. A further study by Snow & Hillier (2019), who used a magnetic field inclined with respect to the direction of propagation, reveal long-lived intermediate (Alfvén) shocks within the slow-mode shock with a shock transition from above to below the Alfvén speed and a reversal of the magnetic field across the shock front. Studying the structure of the shocks requires very good spatial resolution and will be the subject of our future investigations.

Acknowledgments

This work was supported by the Spanish Ministry of Science through the project AYA2014-55078-P and the US National Science Foudation. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. It contributes to the deliverable identified in FP7 European Research Council grant agreement ERC-2017-CoG771310-PI2FA for the project “Partial Ionization: Two-fluid Approach”. The author(s) wish to acknowledge the contribution of Teide High-Performance Computing facilities to the results of this research. TeideHPC facilities are provided by the Instituto Tecnológico y de Energías Renovables (ITER, SA), http://teidehpc.iter.es

References

All Figures

Fig. 1. Parameters of background model atmosphere as function of height. Panel a: gas pressure of neutrals (blue), charges (orange), and magnetic pressures corresponding to the S magnetic field profile (green) and B profile (red). Panel b: number density of neutrals (blue), and charges (orange). For comparison, corresponding number densities from VALC model atmosphere are plotted in green and red. Panel c: temperature. Panel d: horizontal magnetic field, Bx0, of the two background profiles, S and B. Here and below, zero of z axis is considered to be the same as in VALC atmosphere.

In the text

	Fig. 2. Height dependence of characteristic frequencies calculated for background atmospheric model. The blue and orange lines are ion-cyclotron frequencies corresponding to the S and B magnetic field profiles, red and green lines are ion-neutral and neutral-ion collision frequencies, and violet is the highest wave frequency used is our work, corresponding to the wave period of 1 s.
In the text

Fig. 3. Height dependence of oscillations in velocity and temperature of charges and neutrals, and in magnetic field for fixed moment of time, obtained from numerical solution of two fluid equations for magnetic field profile S. Panels a and b: initial wave amplitude factor A = 1; panels c and d are for A = 100. The wave period is P = 20 s for panels a and c, and P = 7.5 s for panels b and d. Orange lines are for the parameters corresponding to charges, and blue lines are those for neutrals. We note that shock overshooting observed in panels c and d are numerical artifacts since we ran the simulation without any shock capturing algorithms in this case.

In the text

	Fig. 4. Height dependence of oscillations in velocity and temperature of charges and neutrals, and in magnetic field for fixed moment in time, for initial wave amplitude factor A = 1. The format of the figure is the same as for Fig. 3. Panels a and b: S magnetic field profile; panels c and d are for the B profile. The wave period is P = 5 s for panels a and c, and P = 1 s for panels b and d.
In the text

	Fig. 5. Decoupling between charges and neutral velocity as function of height for magnetic field profile S. Left panel: wave period 20 s, amplitude factors A = 1 (blue), A = 10 (orange), A = 100 (green). Right panel: wave period 1 s, amplitude factors A = 0.5 (blue), A = 1 (orange), A = 2 (green).
In the text

	Fig. 6. Decoupling between charges and neutral velocity as function of height for magnetic field profile S and wave amplitude factor A = 1. Wave periods P = 20 s (violet), P = 7.5 s (red), P = 5 s (green), P = 2.5 s (orange), and P = 1 s (blue).
In the text

	Fig. 7. Decoupling between charges and neutral velocity as function of height for amplitude factor A = 1 and wave period P = 1 s. The dependence on the magnetic field profiles is shown with different curves: S (blue) and B (orange).
In the text

	Fig. 8. Decoupling between charges and neutral velocity as function of height for S magnetic field profile and amplitude factor A = 1. Blue line: nonlinear simulation; orange line: linear simulation; green line: linear simulation with artificially increased collisional parameter α.
In the text

	Fig. 9. Time average of perturbed variables as function of height for simulation with S magnetic field profile, wave period P = 1 s and amplitude factor A = 2. From top to bottom: temperatures of neutrals (orange) and charges (blue), magnetic field, density of charges, density of neutrals, velocity of neutrals (orange) and charges (blue).
In the text

Fig. 10. Dissipative terms in temperature equations (Eq. (26)) as function of height for simulation with S magnetic field profile, wave period P = 1 s, and amplitude factor A = 2. Top panel: thermal exchange term (TE, orange line) and frictional heating term (FH, blue line) for neutrals. Middle panel: same for charges. Bottom panel: total energy exchange terms for charges (blue), neutrals (orange). We note that apparent oscillations in the bottom panel are numerical artifacts due to averaging the signal with relatively low temporal resolution.

In the text

	Fig. 11. Running average of perturbed temperature (left) and magnetic field (right) as function of time, for fixed height z = 1 Mm for simulation with S magnetic field profile, wave period P = 1 s. The averaging is taken over three consecutive wave periods. The amplitude factor is A = 0.5 (blue), A = 1 (green), and A = 2 (red). Temperature of neutrals and charges closely match together.
In the text

Fig. 12. Solutions of non-dimensional dispersion relation without stratification, Eq. (42). Upper and lower panels: real and imaginary part of k, scaled with $\sqrt{a}/\omega$ as a function of the ratio between the collisional and wave frequencies, scaled with ρ0. Solutions labeled as “S” (dashed lines) are for S magnetic field; solutions labeled “B” (solid lines) are for B magnetic field. The vertical dotted lines mark values of ω = νin (Fc), ω = νni (Fn), and ω = 2π/1 (F1) for the 1 s period wave.

In the text

	Fig. 13. Solutions of dispersion relation (36) for two-fluid linearized equations. Upper and lower panels: real and imaginary part of k as a function of height, calculated for the magnetic field profile S and 5 s period.
In the text

Fig. 14. Comparison between numerical solution in linear regime (black lines) and analytical solution of two-fluid equations (red dotted line). Individual panels show snapshots of the velocity of charges as a function of height at fixed time moments in the stationary regime of the simulations. Panels from left to right, from top to bottom are for wave periods of 1, 5, 7.5, and 20 s, and S magnetic field profile. It is important to note that the analytical solution slightly overshoots the numerical solution, the effect being more pronounced toward the upper part of the domain. This is due to the fact that the analytical solution is an approximate solution.

In the text

Fig. 15. Left panel: imaginary part of k as function of wave period (horizontal axis) and height (vertical axis) obtained after solving dispersion relation for two-fluid equations (Eq. (36)) for S magnetic field profile. Solution #3 from Fig. 13 is used. The period varies between 1 s and 20 s (the minimum and maximum values used in the simulations). The solid, black line represents the contour of kI = 0. Right panel: imaginary part of k as function of background magnetic field at base of atmosphere (horizontal axis) and height (vertical axis). The wave period is fixed to 5 s.

In the text

Fig. 16. Comparison between numerical solution in linear regime (black lines) and nonlinear regime (red dotted line) with A = 10. Individual panels show snapshots of the velocity of charges as a function of height at fixed time moments in the stationary regime of the simulations. Panels from left to right, from top to bottom are for wave periods of 1, 5, 7.5, and 20 s, and S magnetic field profile.

In the text