Impulsively generated two-ﬂuid magnetoacoustic-gravity waves: Solar chromosphere heating and plasma outﬂows

Context. We use a two-ﬂuid model to study the heating of the solar chromosphere by magnetoacoustic and magnetoacoustic-gravity waves. In the model, we include energy dissipation as a result of ion–neutral collisions. Aims. The aim of this paper is to study impulsively generated two-ﬂuid magnetoacoustic and magnetoacoustic-gravity waves and to quantify their contribution to chromosphere heating and the generation of plasma outﬂows. Methods. We consider a 2D model of the gravitationally stratiﬁed and partially ionized solar atmosphere that is permeated by a vertical magnetic ﬁeld. To describe the dynamics of the atmosphere, we use a set of two-ﬂuid equations which we solve numerically with the use of the JOANNA code. Results. We show that large-amplitude impulsively generated magnetoacoustic-gravity waves can e ﬃ ciently heat the chromosphere and generate plasma outﬂows in the low solar corona. The chromosphere is heated by ion–neutral collisions, which are most e ﬀ ective at the top of this atmospheric layer. Wider and larger amplitude pulses heat the atmosphere more e ﬀ ectively and generate faster plasma outﬂows. Conclusions. Large-amplitude, impulsively generated two-ﬂuid magnetoacoustic-gravity waves have the potential to contribute to the solar chromosphere heating and plasma outﬂows in the low corona.


Introduction
The solar atmosphere can be divided into three main plasma layers based on the differences in physical quantities such as temperature, plasma density, and ionization degree. The lowest region, called the photosphere, extends up to 0.5 Mm and has a relatively low temperature of about 5600 K at its bottom which decreases with height to a minimum of about 4300 K at its top. Higher up is the 2 Mm wide chromosphere with a temperature that rises to about 10 4 −10 5 K. The low temperature of the bottom atmospheric layers makes it weakly ionized with an ionization degree of only about 10 −4 at the top of the photosphere, but it increases with height from there onward as the temperature also increases (Khomenko et al. 2014). The corona, which is the top solar atmosphere layer, has an average temperature of 1-3 MK. As a result of this high temperature, the plasma in this region is essentially fully ionized.
The solar atmosphere is penetrated by a diversity of waves. Among them, magnetohydrodynamic (MHD) waves can be distinguished, including fast and slow magnetoacoustic waves (MAWs), Alfvén waves, and the entropy mode (Alfvén 1942;Thomas 1983;Nakariakov & Verwichte 2005;Roberts 2006;Ballester et al. 2018). These waves may form shocks that locally heat the atmosphere (e.g., Ulmschneider et al. 1978;Carlsson & Stein 1995). As the chromosphere is only partially ionized, neutrals can play an important role in heating this layer. One of the different dissipation mechanisms is ion-neutral collisions which can lead to wave damping and the thermalization of wave energy (Erdélyi & James 2004;Khodachenko et al. 2004;Forteza et al. 2007). The presence of neutral particles can thus result in the damping of Alfvén waves in chromospheric spicules (Zaqarashvili et al. 2011). Alfvén waves were also studied in the context of collisional damping (Leake et al. 2005;Soler 2013; Kuźma et al. 2020). These collisional energy dissipation mechanisms constitute an important component in the chromosphere heating.
Internal gravity (henceforth gravity) waves are disturbances in a gravitationally stratified medium. They are anisotropic as they are unable to propagate along the gravity action. Wave periods of these oscillations are higher than the Brunt-Väisälä cut-off period (Moore & Spiegel 1964). The first studies regarding the heating of the solar atmosphere by gravity waves were performed by Whitaker (1963), following the suggestions of Hines (1960) who showed that these waves significantly contribute to the transfer of energy, which should be included in heating models. The energy transfer between gravity and MAWs may be important in the context of atmosphere heating. It was shown that gravity waves can heat the lower chromosphere due to radiative damping (Lighthill 1967). The first piece of observational evidence of internal waves in the solar atmosphere was provided by Schmieder (1976) and Cram (1978). Studies by Mihalas & Toomre (1981, 1982 showed that internal wave dissipation dominates the mid-chromosphere, transferring all energy to these heights. Studies of ion magnetoacoustic-gravity and neutral acoustic-gravity waves by Vigeesh et al. (2017) show local heating of magnetic flux tubes.
In the previous paper, Niedziela et al. (2021, henceforth Paper I) show that triggered photospheric MAWs can heat the chromosphere through ion-neutral collisions. The authors adopted a one-dimensional (1D) model which does not enable internal gravity waves to propagate. As a result, this model requires a refinement by extending it into 2D geometry, which is the goal of the present paper. In this way, internal gravity waves naturally appear in the model in which we focus on chromospheric heating by ion-neutral collisions and related plasma outflows that are generated in the solar corona.
This paper is organized as follows. The two-fluid equations, the equilibrium model, and the initial pulse are presented in Sect. 2. The setup of the numerical simulations is described in Sect. 3 and the results of these simulations are described in Sect. 4. This paper is concluded by a discussion and a summary of the results in Sect. 5.

Two-fluid equations
In our physical model, we consider the gravitationally stratified solar atmosphere. The plasma in the chromosphere is partially ionized and therefore contains both ions and neutrals. We assume that the chromospheric plasmas consist of protons and neutrals, and the ions are made of hydrogen atoms (i.e. they have a charge of +1). A realistic model to describe the plasma dynamics in lower atmosphere layers is based on using the two-fluid equations which can be split into those for neutrals and those for ions + electrons (Zaqarashvili et al. 2011;Leake et al. 2014;Oliver et al. 2016;Maneva et al. 2017;Popescu Braileanu et al. 2019). The set of equations for neutrals (hereafter indicated by the index n) is given by the following: while ions ( i ) and electrons ( e ) are treated as a separate single fluid, which is described by the following equations: Here, V i,n denote the velocities of ions and neutrals, i,n are the mass densities, p ie,n are the gas pressures with p ie = p i + p e = 2p i , E i,n correspond to the total energy densities, I is the identity matrix, B is the magnetic field, and µ is the magnetic permeability. Constant quantities such as the specific heat ratio, γ, and gravity vector, g = [0, −g, 0], have the magnitudes of γ = 5/3 and g = 274.78 m s −2 . The interaction between the two plasma species takes place through ion-neutral collisions and it depends on the ion-neutral friction coefficient, given as (Braginskii 1965;Oliver et al. 2016) Here, k B is the Boltzmann constant, σ in = 1.4 × 10 −15 cm 2 corresponds to the quantum collision cross section taken from the model of Vranjes & Krstic (2013), and m i and m n are the atomic masses of ions and neutrals equal to the proton mass, respectively. This interaction results in additional source terms in the equations such as the energy exchange terms Q i,n , which are defined as (Meier & Shumlak 2012;Oliver et al. 2016) For simplicity reasons, we neglected all nonideal and nonadiabatic terms.

Magnetohydrostatic equilibrium
It is assumed that the atmosphere is initially (at t = 0 s) in hydrostatic equilibrium (V i,n = 0) in which state the pressure gradient and the gravity forces balance each other: Hence, background gas pressures and mass densities can be derived from (e.g., Murawski et al. 2015) i,n (y) = where p i0 = 10 −5 Pa and p n0 = 3 × 10 −7 Pa are the gas pressures at the reference level, taken at y = y r = 50 Mm, and denote the ion and neutral-pressure scale heights, respectively. Figure 1 (top left) shows vertical variation of the equilibrium temperature (Avrett & Loeser 2008). The minimum temperature of T = 4341 K takes place in the low chromosphere, at y ≈ 0.6 Mm. In the middle and the upper chromosphere, T (y) slightly increases, but the fast growth of the initial temperature takes place at the transition region where the temperature rises from about 7 × 10 3 K to 3 × 10 5 K at y = 3 Mm. The corona is the hottest region reaching T ≈ 1 MK at y = 20 Mm (not shown). The equilibrium profiles of the ion gas pressures and ion mass densities, given by Eqs. (14) and (15), are overlaid by a vertical magnetic field B = B y e y of magnitude B y = 10 G, unless stated otherwise. Henceforth, it is assumed that the initial temperatures for both species are equal, T i = T n = T . A32, page 2 of 11

Bulk cut-off periods and plasma β
In the solar atmosphere, internal-gravity waves can propagate only if their wave periods are larger than the bulk Brunt-Väisälä period, P BV , which is given as (e.g., Vigeesh et al. 2021) Here, is the bulk-pressure scale height, and is the bulk mass-density scale height. Generally, P BV attains values larger than about 100 s ( Fig. 1, top right), which are much larger than the ion-neutral collision timescale in the chromosphere (Khomenko et al. 2014). Acoustic waves can propagate freely in the medium if their wave period P is smaller than the total acoustic cut-off period, given as (Deubner & Gough 1984) P ac = 4πΛ Here, C S = γ(p i + p n )/( i + n ) is the bulk sound speed. In the opposite case, mainly for P > P ac , acoustic waves are evanescent (Lamb 1932).
We define now the bulk plasma β as the ratio of total thermal pressure to magnetic pressure and it is expressed as The total plasma β is an important coefficient that provides us information about the impact of the magnetic and thermal pressures in the solar atmosphere. We note that β 1 in the photosphere and bottom chromosphere (Fig. 1, bottom right), which implies that the total thermal pressure dominates over its magnetic counterpart, and slow and fast MAWs are strongly coupled and indistinguishable. The magnitude of the plasma β decreases with height, resulting in a weaker coupling of the MAWs. The solar corona is the region of plasma β < 1 from which we infer that slow and fast MAWs are distinguishable and the coronal plasma is governed by the Lorentz force, while the gas pressure force is less important (Nakariakov & Verwichte 2005).

Initial pulse
To perturb magnetohydrostatic equilibrium, we used the Gaussian pulse in the vertical components of the ion and neutral velocities, in other words where A is the amplitude of the pulse, w corresponds to its width, and (0, y 0 ) is the point from which the pulse is launched. We

Numerical simulations
We performed numerical simulations with the use of the JOANNA code which solves the initial-boundary-value problem for the two-fluid equations ). The 2D numerical box is specified as −2.56 Mm ≤ x ≤ 2.56 Mm along the horizontal (x−) direction with the cell size ∆x = 2.5 km. The region along the vertical (y−) direction is covered at the bottom by a fine uniform grid zone with cells of size ∆y = 2.5 km. In the case of y 0 = 0 Mm, the fine grid zone starts at y = −2 Mm and ends at y = 3.12 Mm. For y 0 = 0.5 Mm, the fine grid zone is limited by min (

Numerical tests
The results of each numerical simulation are subject to certain numerical errors. We, therefore, started the simulations by study- ing the role of the background noise produced by the JOANNA code which helped us to verify the numerical accuracy of the simulation setup. The top panel of Fig. 2 displays time-distance plots for the ratio of the vertical component of the ion velocity and the bulk sound speed. We note that max |V iy /C S | ≈ 4 × 10 −2 and this maximum occurs in the corona, highlighting that the numerical error is not dynamically relevant; from this, we deduced that V iy is small in comparison to the local sound speed in this region Figure 2 (bottom) shows δT i /T 0 = (T i − T 0 )/T 0 from which we conclude that max |δT i /T 0 | ≈ 4 × 10 −3 occurs in the top chromosphere and it is very small. In addition, these values also decay with time, achieving 2 × 10 −5 of their initial values after 1750 s. Hence, we infer that the numerical errors seen in both panels are negligibly small, and the chosen numerical grid and discretization scheme do not essentially affect the numerical results.

A pulse launched from the bottom of the photosphere
First, we consider a pulse with amplitude A = 2 km s −1 . Figure 3 shows time-distance plots for V iy (top) and δT i /T 0 (bottom), collected along the y axis at x = 0. Ions reach a maximum velocity of V iy ≈ 30 km s −1 . Steep wave profiles are formed in the chromosphere and result in plasma heating through ion-neutral collisions. Hence, we observe the spreading of the initial pulse into a train of waves. At t ≈ 300 s, a part of the signal is reflected from the transition region. Thus, the upward propagating waves travel with an average speed of about 7.14 km s −1 . This value is comparable to the bulk sound speed at the bottom of the A32, page 4 of 11 chromosphere, C S ≈ 6.9 km s −1 at y = 0.7 Mm, and to the bulk Alfvén speed at the middle of this layer, V A = B y / µ( i + n ) ≈ 7 km s −1 at y = 1.4 Mm. A reduction in the temperature of the chromosphere is observed during the initial phase (bottom) due to a rarefaction. After t = 300 s, the temperature of the chromosphere essentially grows. The pulse was initially (at t = 0 s) launched from the bottom of the photosphere (Fig. 4, top left). As the signal spreads in space, its amplitude decreases to A = 1.74 km s −1 at t = 50 s (Fig. 4, top right). At t = 100 s, the signal, which continues to move upward, increases in amplitude (Fig. 4, bottom left) due to the decreasing background density. At t = 200 s, the signal already passed the transition region, resulting in a significant wave steepening (Fig. 4, bottom right).

Small amplitude pulse
In this subsection, we consider the case of y 0 = 0.5 Mm and A = 2 km s −1 . Figure 5 presents time-distance plots for V iy collected at x = 0 Mm (left) and x = 0.5 Mm (right). The left panel corresponds to MAWs, while the right panel is associated with MAWs and gravity waves. Plots of V iy are drawn up to y = 5 Mm, while δT i /T 0 and V iy − V ny up to 2 Mm to show more details, which would hardly be discernible if these plots were redrawn for the whole domain. In both panels, the signals hit the transition region at t ≈ 200 s and experience a partial reflection from there. The vertical component of the A32, page 5 of 11 A&A 668, A32 (2022) ion velocity, V iy , collected for x = 0 Mm reaches a higher value of max(V iy ) = 17 km s −1 , than for x = 0.5 Mm, at which max(V iy ) = 10 km s −1 . Shocks that are formed at the transition region heat the atmosphere by ion-neutral collisions. The heating of the atmosphere by MAWs is more pronounced with larger maximum δT i /T 0 ≈ 0.7 values compared to magnetoacousticgravity waves (MAGWs). Figure 5 (bottom) shows log of the absolute value of velocity drift, log |δV| = log |V iy − V ny |, which grows with height and reaches its maximum at the top of the chromosphere. It shows that the top of the chromosphere is essentially heated by collisions (see Eqs. (11) and (12)). The highest value of log |δV| occurs during the initial phase when the pulse reaches the transition region at t = 200 s. This is in agreement with Fig. 5, which illustrates that the plasma is most strongly heated at t = 200 s. The minimum value of log |δV| takes place below y = 1 Mm from which it follows that ions and neutrals are strongly coupled and propagate with almost the same speed. As a result of the ion-neutral collisions, the wave energy is dissipated and converted into heat. This dissipation mechanism is most effective for the largest values of log |δV| and it A32, page 6 of 11 can compensate the radiative and thermal losses in the chromosphere. The results obtained in Paper I also show that the highest values of log |δV| take place at the top of the chromosphere. Figure 6 illustrates V iy (top) and δT i /T 0 (bottom), averaged over time, that is specified as where f = V iy or δT i /T 0 , t 1 = 500 s, and t 2 = 2000 s. For MAWs (left), V iy starts to increase from the bottom of the chromosphere y = 0.5 Mm and it reaches its maximum of about 0.04 km s −1 at y ≈ 0.7 Mm. A minimum of V iy = −0.25 km s −1 is observed at y = 5 Mm. On the contrary, the panel for MAGWs (right) shows an increase in V iy from the middle of the chromosphere y = 1.1 Mm up to the transition region. For MAWs, upflows (V iy > 0) are observed within the range of 0.5 Mm < y < 1.5 Mm, while for MAGWs they are observed above y ≈ 1.5 Mm. Higher up in the solar corona, such upflows may result in the origin of solar wind. These results are in agreement with Paper I. It is apparent that MAWs heat the chromosphere more significantly than MAGWs. As a matter of fact, for MAWs, δT i /T 0 reaches a value of about 0.4 at the top of the chromosphere. For the left panel (MAWs), a minimum of δT i /T 0 is observed at y 0 = 0.5 Mm, which corresponds to the launching height in this case. On the other hand, in the RHS panel (MAGWs), we observe temperature oscillations at the launching level, while a temperature minimum is present at y 0 = 1 Mm. For both types of waves, the top of the chromosphere is heated the most, which correlates with δV iy (Fig. 5, bottom). We note that the temperature in the photosphere is reduced since no waves are seen there (Fig. 5, bottom). A comparison with the results from Paper I shows that the average heating in the 2D model is significantly smaller than that obtained with the 1D model. Figure 7 shows the kinetic energy flux. This quantity is given by the equation The top panels clearly show that the greatest energy loss takes place in the chromosphere with the highest value of max(F Ecsi ) = 10 5 erg cm −2 s −1 for the LHS panel and max(F Ecsi ) = 10 6 erg cm −2 s −1 for the RHS panel. Hence, we observe higher energy losses for MAWs than MAGWs. The bottom panels show F Ecsi averaged over time. We note local minima occurring in the middle of the chromosphere for both types of waves. Above the transition region, we observe a plateau of F Ecsi . The obtained kinetic energy flux values for MAWs are too small to fulfill radiative losses in the chromosphere by two orders of magnitude (Withbroe & Noyes 1977). For MAGWs, it is a difference of three orders of magnitude. A comparison of F Ecsi in the corona similarly shows that these are not sufficient values to match the energy losses. These differences may be the result of the absence of nonlinear Alfvén waves.

Fourier power spectra
Figure 8 displays wave periods from the Fourier power spectrum V iy (x = 0, y, t) (left) and V iy (x = 0.5, y, t) (right). The wave cutoffs that take place in stratified media can be used to determine if the range of periods corresponds to propagating or evanescent waves. Observational variations of the cutoff in the solar atmosphere were studied by Wiśniewska et al. (2016) and Kayshap et al. (2018). For x = 0 Mm, the main period P ≈ 200 s. The wave periods are only higher than P AC in the lower photosphere, since max(P AC ) = 220 s at y = 0.25 Mm, and in the transition region within the range of 2.1 Mm < y < 2.2 Mm. Hence, the waves are evanescent outside of these regions. For x = 0.5 Mm, the dominant wave period is close to P ≈ 300 s within the range of 0 Mm < y < 0.6 Mm. Higher up, for 0.6 Mm < y < 1.5 Mm, the signal propagates with P ≈ 230 s. At y = 1.5 Mm, the main period falls off to P ≈ 150 s. Hence, the waves can propagate throughout the photosphere and the low chromosphere. A comparison of these wave periods with the observational data of Wiśniewska et al. (2016) and Kayshap et al. (2018) reveals an agreement at some points. Hence, the similarity between the numerical results and the observational data confirms that ionneutral collisions are an efficient mechanism of energy release and our results can be used to determine the structure of the solar atmosphere.

Parametric studies
We performed some parametric studies to investigate how the pulse width affects plasma heating and the generation of outflows. Decreasing the pulse width from w = 0.3 Mm to w = 0.2 Mm results in smaller V iy at a given height (Fig. 9, top). For a narrower pulse width, max(δT i /T 0 ) falls off by a factor of two (bottom). Thus, the chromosphere is not heated by narrower pulses as much. Figure 10 illustrates the maximum of V iy (x = 0, y, t) for different pulse parameters. A higher absolute amplitude, |A| (top), and a higher width of the pulse (bottom) both result in an increase of max(V iy ), which is in agreement with the 1D results of Paper I. However, the 1D model exhibited a linear correlation between the width of the pulse and the maximum outflow velocity, while the increase is not linear in the 2D model.
As a result of increasing the absolute value of the amplitude of the initial pulse from A = 2 km s −1 to A = 5 km s −1 , the signal in V iy , which was collected at a given height, is more pronounced for larger amplitudes (Fig. 11, left). The maximum values of δT i /T 0 also increased to about 1.6 (right) compared to the results for smaller amplitude max(δT i /T 0 ) = 0.7 (Fig. 5, middle left) at the top of the chromosphere. However, comparison with the results for y 0 = 0 Mm (Fig. 3) shows that even though the amplitude is larger, the signal in V iy is smaller. This is a result of the lower y 0 value which affects the pressure-scale height. A32, page 8 of 11   9. Time-distance plots for V iy (top) and δT i /T 0 (bottom), collected along y for x = 0 Mm, for the pulse of its amplitude A = 2 km s −1 , and width w = 0.2 Mm, launched from y = y 0 = 0.5 Mm. Figure 12 shows parametric studies for different background magnetic field magnitudes. A magnetic field of magnitude B y = 5 G (top-left) results in max(δT i /T 0 ) ≈ 0.5. We note that plasma is heated and cooled within the whole atmosphere. Increasing the magnitude of the magnetic field to B y = 15 G (top right) results in higher max(δT i /T 0 ) ≈ 1.2. In this case, after t = 1200 s, no more cooling is observed at the top of the chromosphere. Magnetic fields of even higher magnitudes, such as B y = 20 G (bottom left) and B y = 25 G (bottom right), cause an increase in max(δT i /T 0 ) to 2.1 and 2.8, respectively. After t = 1200 s, the top of the chromosphere is hardly cooled at all. A32, page 9 of 11 A&A 668, A32 (2022) Fig. 11. Time-distance plots for V iy (left) and δT i /T (right), collected along y for x = 0 Mm, for the pulse of its amplitude A = 5 km s −1 , and width w = 0.3 Mm, launched from y = y 0 = 0.5 Mm.  Figure 13 presents the relative perturbed ion temperature averaged over time and height, H, versus B y . We define this quantity as where y 0 = 0.5 Mm and y 1 = 1.9 Mm. It is clearly seen that higher magnitudes of the vertical magnetic field result in more heating of the photosphere and chromosphere.

Summary and conclusion
We performed 2D numerical simulations of impulsively generated neutral acoustic-gravity and ion magnetoacoustic-gravity two-fluid waves in the partially ionized lower solar atmosphere. These waves are triggered by the initially launched pulse in the vertical components of both the ion and neutral velocities, and they heat the upper chromosphere as a result of ion-neutral collisions. Increasing the launching height results in a smaller temperature increase in the chromosphere. This results from the A32, page 10 of 11 pressure scale height according to which the amplitude increases e times with height. Moreover, max(V iy ) and max(δT i /T 0 ) are higher for vertically propagating magnetoacoustic waves than for magnetoacoustic-gravity waves which propagate obliquely to the gravity action, that is to say horizontally. The heating of the chromosphere is strongly correlated with the velocity drift which attains its largest values at the top of this region. The same results were obtained in Paper I where the top chromosphere was heated as a result of ion-neutral collisions. Similarly to the plasma heating, the plasma outflow is also larger for magnetoacoustic waves than for magnetoacoustic-gravity waves. Parametric studies show that a narrower pulse results in less chromosphere heating and lower plasma outflows, which results from less energy being released by a narrower pulse. On the other hand, increasing the absolute value of the amplitude from A = 2 km s −1 to A = 5 km s −1 leads to more significant heating and a more substantial generation of plasma outflows. The parametric studies in Paper I for the 1D model have also shown that a narrower pulse leads to less heating of the atmosphere and smaller plasma outflows and that increasing the amplitude leads to higher V iy and δT i /T 0 . As a final conclusion, we note that larger-amplitude magnetoacoustic waves excited by a single pulse result in stronger chromospheric heating and the generation of plasma outflows for both 1D and 2D models.