© R. Niedziela et al. 2022

1. 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⁴ − 10⁵ K. The low temperature of the bottom atmospheric layers makes it weakly ionized with an ionization degree of only about 10⁻⁴ 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.

2. Numerical model of the solar atmosphere

2.1. 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:

$$\begin{array}{cc}& \frac{\partial {\varrho }_{\mathrm{n}}}{\partial t}+\mathrm{\nabla }·({\varrho }_{\mathrm{n}}{\mathit{V}}_{\mathrm{n}})=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\partial ({\varrho }_{\mathrm{n}}{\mathit{V}}_{\mathrm{n}})}{\partial t}+\mathrm{\nabla }·({\varrho }_{\mathrm{n}}{\mathit{V}}_{\mathrm{n}}{\mathit{V}}_{\mathrm{n}}+p_{\mathrm{n}}\mathit{I})={\varrho }_{\mathrm{n}}\mathit{g}-{\alpha }_{\mathrm{c}}({\mathit{V}}_{\mathrm{n}}-{\mathit{V}}_{\mathrm{i}}),\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial E_{\mathrm{n}}}{\partial t}+\mathrm{\nabla }·[(E_{\mathrm{n}}+p_{\mathrm{n}}){\mathit{V}}_{\mathrm{n}}]={\varrho }_{\mathrm{n}}\mathit{g}·{\mathit{V}}_{\mathrm{n}}+Q_{\mathrm{n}},\hfill \end{array}$$(3)

$$\begin{array}{cc}& E_{\mathrm{n}}=\frac{{\varrho }_{\mathrm{n}}{\mathit{V}}_{\mathrm{n}}^2}{2}+\frac{p_{\mathrm{n}}}{\gamma -1},\hfill \end{array}$$(4)

while ions (_i) and electrons (_e) are treated as a separate single fluid, which is described by the following equations:

$$\begin{array}{cc}& \frac{\partial {\varrho }_{\mathrm{i}}}{\partial t}+\mathrm{\nabla }·({\varrho }_{\mathrm{i}}{\mathit{V}}_{\mathrm{i}})=0,\hfill \end{array}$$(5)

$$\begin{array}{cc}& \frac{\partial ({\varrho }_{\mathrm{i}}{\mathit{V}}_{\mathrm{i}})}{\partial t}+\mathrm{\nabla }·({\varrho }_{\mathrm{i}}{\mathit{V}}_{\mathrm{i}}{\mathit{V}}_{\mathrm{i}}+p_{\mathrm{ie}}\mathit{I})={\varrho }_{\mathrm{i}}\mathit{g}+\frac{1}{\mu }(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}+{\alpha }_{\mathrm{c}}({\mathit{V}}_{\mathrm{n}}-{\mathit{V}}_{\mathrm{i}}),\hfill \end{array}$$(6)

$$\begin{array}{cc}& \frac{\partial E_{\mathrm{i}}}{\partial t}+\mathrm{\nabla }·[(E_{\mathrm{i}}+p_{\mathrm{ie}}+\frac{{\mathit{B}}^2}{2\mu }){\mathit{V}}_{\mathrm{i}}-\frac{\mathit{B}}{\mu }({\mathit{V}}_{\mathrm{i}}·\mathit{B})]={\varrho }_{\mathrm{i}}\mathit{g}·{\mathit{V}}_{\mathrm{i}}+Q_{\mathrm{i}},\hfill \end{array}$$(7)

$$\begin{array}{cc}& E_{\mathrm{i}}=\frac{{\varrho }_{\mathrm{i}}{\mathit{V}}_{\mathrm{i}}^2}{2}+\frac{p_{\mathrm{ie}}}{\gamma -1}+\frac{{\mathit{B}}^2}{2\mu },\hfill \end{array}$$(8)

$$\begin{array}{cc}& \frac{\partial \mathit{B}}{\partial t}=\mathrm{\nabla }\times ({\mathit{V}}_{\mathrm{i}}\times \mathit{B}),\mathrm{\nabla }·\mathit{B}=0.\hfill \end{array}$$(9)

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⁻². 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)

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{c}}=\frac{4}{3}\frac{{\sigma }_{\mathrm{in}}}{m_{\mathrm{i}}+m_{\mathrm{n}}}\sqrt{\frac{8k_{\mathrm{B}}}{\pi }(\frac{T_{\mathrm{i}}}{m_{\mathrm{i}}}+\frac{T_{\mathrm{n}}}{m_{\mathrm{n}}})}{\varrho }_{\mathrm{n}}{\varrho }_{\mathrm{i}}.\end{array}$$(10)

Here, k_B is the Boltzmann constant, σ_in = 1.4 × 10⁻¹⁵ cm² 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)

$$\begin{array}{cc}& Q_{\mathrm{i}}={\alpha }_{\mathrm{c}}[\frac{1}{2}{({\mathit{V}}_{\mathrm{i}}-{\mathit{V}}_{\mathrm{n}})}^2+3\frac{k_{\mathrm{B}}}{m_{\mathrm{i}}+m_{\mathrm{n}}}(T_{\mathrm{n}}-T_{\mathrm{i}})],\hfill \end{array}$$(11)

$$\begin{array}{cc}& Q_{\mathrm{n}}={\alpha }_{\mathrm{c}}[\frac{1}{2}{({\mathit{V}}_{\mathrm{i}}-{\mathit{V}}_{\mathrm{n}})}^2+3\frac{k_{\mathrm{B}}}{m_{\mathrm{i}}+m_{\mathrm{n}}}(T_{\mathrm{i}}-T_{\mathrm{n}})].\hfill \end{array}$$(12)

For simplicity reasons, we neglected all nonideal and nonadiabatic terms.

2.2. Initial conditions

2.2.1. 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:

$$\begin{array}{c}\hfill -\mathrm{\nabla }{p}_{\mathrm{i},\mathrm{n}}+{\varrho }_{\mathrm{i},\mathrm{n}}\mathit{g}=\mathbf{0}.\end{array}$$(13)

Hence, background gas pressures and mass densities can be derived from (e.g., Murawski et al. 2015)

$$\begin{array}{cc}& p_{\mathrm{i},\mathrm{n}}(y)=p_{\mathrm{i},\mathrm{n}0}exp(-{\displaystyle {\int }_{y_r}^y}\frac{\mathrm{d}y}{{\mathrm{\Lambda }}_{\mathrm{i},\mathrm{n}}(y)}),\hfill \end{array}$$(14)

$$\begin{array}{cc}& {\varrho }_{\mathrm{i},\mathrm{n}}(y)=\frac{p_{\mathrm{i},\mathrm{n}}(y)}{g{\mathrm{\Lambda }}_{\mathrm{i},\mathrm{n}}(y)},\hfill \end{array}$$(15)

where p_i0 = 10⁻⁵ Pa and p_n0 = 3 × 10⁻⁷ Pa are the gas pressures at the reference level, taken at y = y_r = 50 Mm, and

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{i},\mathrm{n}}=\frac{k_{\mathrm{B}}{T}_{\mathrm{i},\mathrm{n}}}{m_{\mathrm{i},\mathrm{n}}g}\end{array}$$(16)

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³ K to 3 × 10⁵ 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_ye_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.

Fig. 1. Variation with height of the equilibrium temperature (top left), the bulk Brunt–Väisälä period (top right), the acoustic cut-off period (bottom left), and the plasma β (bottom right) in the initial equilibrium configuration.

2.3. 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)

$$\begin{array}{c}\hfill P_{\mathrm{BV}}=2\pi \sqrt{\frac{1}{g}(\gamma \mathrm{\Lambda }-H)}.\end{array}$$(17)

Here,

$$\begin{array}{c}\hfill \mathrm{\Lambda }=\frac{k_{\mathrm{B}}(T_{\mathrm{i}}+T_{\mathrm{n}})}{(m_{\mathrm{i}}+m_{\mathrm{n}})g}\end{array}$$(18)

is the bulk-pressure scale height, and

$$\begin{array}{c}\hfill H=\left|\frac{{\varrho }_{\mathrm{i}}+{\varrho }_{\mathrm{n}}}{{({\varrho }_{\mathrm{i}}+{\varrho }_{\mathrm{n}})}_{,y}}\right|\end{array}$$(19)

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)

$$\begin{array}{c}\hfill P_{\mathrm{ac}}=\frac{4\pi \mathrm{\Lambda }}{C_{\mathrm{S}}\sqrt{1+2\frac{\mathrm{d}\mathrm{\Lambda }}{\mathrm{d}y}}}·\end{array}$$(20)

Here, $C_{\mathrm{S}}=\sqrt{\gamma (p_{\mathrm{i}}+p_{\mathrm{n}})/({\varrho }_{\mathrm{i}}+{\varrho }_{\mathrm{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

$$\begin{array}{c}\hfill \beta =\frac{p_{\mathrm{ie}}+p_{\mathrm{n}}}{B^2/2\mu }·\end{array}$$(21)

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).

2.3.1. Initial pulse

To perturb magnetohydrostatic equilibrium, we used the Gaussian pulse in the vertical components of the ion and neutral velocities, in other words

$$\begin{array}{c}\hfill V_{\mathrm{iy}}(x,y,t=0)=V_{\mathrm{ny}}(x,y,t=0)=Aexp(-\frac{x^2+{(y-y_0)}^2}{w^2}),\end{array}$$(22)

where A is the amplitude of the pulse, w corresponds to its width, and (0, y₀) is the point from which the pulse is launched. We consider several values of y₀, w, and A, mainly y₀ = 0 Mm and y₀ = 0.5 Mm, w = 0.3 Mm and w = 0.2 Mm, and A = 2 km s⁻¹ and A = −5 km s⁻¹.

3. 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 (Wójcik et al. 2020). 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 Mm, the fine grid zone starts at y = −2 Mm and ends at y = 3.12 Mm. For y₀ = 0.5 Mm, the fine grid zone is limited by min(y) = 0 Mm and max(y) = 5.12 Mm. Above this fine grid zone, a stretched grid is implemented which consists of 32 cells up to y = 40 Mm. Such a small number of cells in the stretched grid zone is required to absorb the incoming signal and essentially to remove reflections from the top boundary at which – along with the bottom boundary – all fluid quantities are set to their equilibrium values. On the left and right sides of the numerical box, all fluid variables were copied from the nearby physical cells into the boundary cells, consisting of open boundary conditions.

4. Numerical results

4.1. Numerical tests

The results of each numerical simulation are subject to certain numerical errors. We, therefore, started the simulations by studying 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⁻² 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

Fig. 2. Time–distance plots for Viy/CS (top) and δTi/T0 (bottom), without any pulse (A = 0) being launched initially (at t = 0 s).

Figure 2 (bottom) shows δT_i/T₀ = (T_i − T₀)/T₀ from which we conclude that max |δT_i/T₀|≈4 × 10⁻³ occurs in the top chromosphere and it is very small. In addition, these values also decay with time, achieving 2 × 10⁻⁵ 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.

4.2. A pulse launched from the bottom of the photosphere

First, we consider a pulse with amplitude A = 2 km s⁻¹. Figure 3 shows time–distance plots for V_iy (top) and δT_i/T₀ (bottom), collected along the y axis at x = 0. Ions reach a maximum velocity of V_iy ≈ 30 km s⁻¹. 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⁻¹. This value is comparable to the bulk sound speed at the bottom of the chromosphere, C_S ≈ 6.9 km s⁻¹ at y = 0.7 Mm, and to the bulk Alfvén speed at the middle of this layer, $V_{\mathrm{A}}=B_{\mathit{y}}/\sqrt{\mu ({\varrho }_{\mathrm{i}}+{\varrho }_{\mathrm{n}})}\approx 7$ km s⁻¹ 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⁻¹ 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).

Fig. 3. Time–distance plots for Viy (top) and δTi/T0 (bottom), collected along y at x = 0 Mm for A = 2 km s−1, w = 0.3 Mm, and y0 = 0 Mm.

Fig. 4. Spatial profiles of Viy expressed in units of km s−1 at t = 0 s (top left), t = 50 s (top right), t = 100 s (bottom left), and t = 200 s (bottom right) for A = 2 km s−1, w = 0.3 Mm, and y0 = 0 Mm.

4.3. A pulse launched from the bottom of the chromosphere

4.4. Small amplitude pulse

In this subsection, we consider the case of y₀ = 0.5 Mm and A = 2 km s⁻¹. 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₀ 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 ion velocity, V_iy, collected for x = 0 Mm reaches a higher value of max(V_iy) = 17 km s⁻¹, than for x = 0.5 Mm, at which max(V_iy) = 10 km s⁻¹. 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.7 values compared to magnetoacoustic-gravity 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 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.

Fig. 5. Time–distance plots for Viy (top), δTi/T0 (middle), and log|Viy−Vny| (bottom), collected along y at x = 0 (left) and x = 0.5 (right), for the pulse of its amplitude A = 2 km s−1, and width w = 0.3 Mm, launched from y = y0 = 0.5 Mm.

Figure 6 illustrates V_iy (top) and δT_i/T₀ (bottom), averaged over time, that is specified as

Fig. 6. Averaged over time ⟨Viy⟩ (top) and ⟨δTi/T0⟩ (bottom), collected along y at x = 0 Mm (left) and x = 0.5 Mm (right), for the pulse parameters corresponding to Fig. 5.

$$\begin{array}{c}\hfill ⟨f⟩=\frac{1}{t_2-t_1}{\displaystyle {\int }_{t_1}^{t_2}}f\mathrm{d}t,\end{array}$$(23)

where f = V_iy or δT_i/T₀, t₁ = 500 s, and t₂ = 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⁻¹ at y ≈ 0.7 Mm. A minimum of V_iy = −0.25 km s⁻¹ 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₀⟩ reaches a value of about 0.4 at the top of the chromosphere. For the left panel (MAWs), a minimum of ⟨δT_i/T₀⟩ is observed at y₀ = 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₀ = 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

$$\begin{array}{c}\hfill F_{\mathrm{Ecsi}}=0.5\varrho V^2{C}_{\mathrm{S}}.\end{array}$$(24)

Fig. 7. Time–distance plots for FEcsi (top) and vertical profiles of ⟨FEcsi⟩ averaged over time (bottom), collected along y at x = 0 Mm (left) and y at x = 0.5 Mm (right), for the pulse parameters corresponding to Fig. 5.

The top panels clearly show that the greatest energy loss takes place in the chromosphere with the highest value of max(F_Ecsi) = 10⁵ erg cm⁻² s⁻¹ for the LHS panel and max(F_Ecsi) = 10⁶ erg cm⁻² s⁻¹ 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.

4.5. 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 ion–neutral collisions are an efficient mechanism of energy release and our results can be used to determine the structure of the solar atmosphere.

Fig. 8. Fourier period P for Viy vs. height, collected along y at x = 0 Mm (left) and x = 0.5 Mm (right), for the pulse of its amplitude A = 2 km s−1, width w = 0.3 Mm, and launched from y = y0 = 0.5 Mm. The observational data of Wiśniewska et al. (2016) and Kayshap et al. (2018) are represented by diamonds and dots with the vertical error bars, respectively.

4.6. 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₀) falls off by a factor of two (bottom). Thus, the chromosphere is not heated by narrower pulses as much.

Fig. 9. Time–distance plots for Viy (top) and δTi/T0 (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 = y0 = 0.5 Mm.

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.

Fig. 10. Maximum value of Viy(x = 0, y, t) vs. the following: |A| for y0 = 0.5 Mm and w = 0.3 Mm (top); and w for y0 = 0.5 Mm and A = 2 km s−1 (bottom).

As a result of increasing the absolute value of the amplitude of the initial pulse from A = 2 km s⁻¹ to A = 5 km s⁻¹, 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₀ also increased to about 1.6 (right) compared to the results for smaller amplitude max(δT_i/T₀) = 0.7 (Fig. 5, middle left) at the top of the chromosphere. However, comparison with the results for y₀ = 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₀ value which affects the pressure-scale height.

Fig. 11. Time–distance plots for Viy (left) and δTi/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 = y0 = 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.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₀) ≈ 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₀) to 2.1 and 2.8, respectively. After t = 1200 s, the top of the chromosphere is hardly cooled at all.

Fig. 12. Time–distance plots for δTi/T, collected along y for x = 0 Mm, for the pulse of its amplitude A = 2 km s−1, and width w = 0.3 Mm, launched from y = y0 = 0.5 Mm, in the case of By = 5 G (top left), By = 15 G (top right), By = 20 G (bottom left), and By = 25 G (bottom right).

Figure 13 presents the relative perturbed ion temperature averaged over time and height, H, versus B_y. We define this quantity as

$$\begin{array}{c}\hfill H=\frac{1}{y_1-y_0}{\displaystyle {\int }_{y_0}^{y_1}}\langle \frac{\delta T_i}{T_0}\rangle \mathrm{d}y,\end{array}$$(25)

Fig. 13. Averaged over time and height relative to perturbed temperature of ions, H vs. By, for the pulse of its amplitude A = 2 km s−1 and width w = 0.3 Mm, launched from y = y0 = 0.5 Mm.

where y₀ = 0.5 Mm and y₁ = 1.9 Mm. It is clearly seen that higher magnitudes of the vertical magnetic field result in more heating of the photosphere and chromosphere.

5. 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 pressure scale height according to which the amplitude increases e times with height. Moreover, max(V_iy) and max(δT_i/T₀) 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⁻¹ to A = 5 km s⁻¹ 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₀.

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.

Acknowledgments

The JOANNA code was developed by Darek Wójcik with some contribution from Luis Kadowaki and Piotr Wołoszkiewicz. This work was done within the framework of the projects from the Polish Science Center (NCN) Grant No. 2020/37/B/ST9/00184. We visualize the simulation data using the VisIt software package (Childs et al. 2012). T.V.Z. was supported by Shota Rustaveli National Science Foundation of Georgia (project FR-21-467). SP acknowledges support from the projects C14/19/089 (C1 project Internal Funds KU Leuven), G.0D07.19N (FWO-Vlaanderen), SIDC Data Exploitation (ESA Prodex-12), and Belspo project B2/191/P1/SWiM.

References

All Figures

	Fig. 1. Variation with height of the equilibrium temperature (top left), the bulk Brunt–Väisälä period (top right), the acoustic cut-off period (bottom left), and the plasma β (bottom right) in the initial equilibrium configuration.
In the text

	Fig. 2. Time–distance plots for Viy/CS (top) and δTi/T0 (bottom), without any pulse (A = 0) being launched initially (at t = 0 s).
In the text

	Fig. 3. Time–distance plots for Viy (top) and δTi/T0 (bottom), collected along y at x = 0 Mm for A = 2 km s−1, w = 0.3 Mm, and y0 = 0 Mm.
In the text

	Fig. 4. Spatial profiles of Viy expressed in units of km s−1 at t = 0 s (top left), t = 50 s (top right), t = 100 s (bottom left), and t = 200 s (bottom right) for A = 2 km s−1, w = 0.3 Mm, and y0 = 0 Mm.
In the text

	Fig. 5. Time–distance plots for Viy (top), δTi/T0 (middle), and log|Viy−Vny| (bottom), collected along y at x = 0 (left) and x = 0.5 (right), for the pulse of its amplitude A = 2 km s−1, and width w = 0.3 Mm, launched from y = y0 = 0.5 Mm.
In the text

	Fig. 6. Averaged over time ⟨Viy⟩ (top) and ⟨δTi/T0⟩ (bottom), collected along y at x = 0 Mm (left) and x = 0.5 Mm (right), for the pulse parameters corresponding to Fig. 5.
In the text

	Fig. 7. Time–distance plots for FEcsi (top) and vertical profiles of ⟨FEcsi⟩ averaged over time (bottom), collected along y at x = 0 Mm (left) and y at x = 0.5 Mm (right), for the pulse parameters corresponding to Fig. 5.
In the text

	Fig. 8. Fourier period P for Viy vs. height, collected along y at x = 0 Mm (left) and x = 0.5 Mm (right), for the pulse of its amplitude A = 2 km s−1, width w = 0.3 Mm, and launched from y = y0 = 0.5 Mm. The observational data of Wiśniewska et al. (2016) and Kayshap et al. (2018) are represented by diamonds and dots with the vertical error bars, respectively.
In the text

	Fig. 9. Time–distance plots for Viy (top) and δTi/T0 (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 = y0 = 0.5 Mm.
In the text

	Fig. 10. Maximum value of Viy(x = 0, y, t) vs. the following: |A| for y0 = 0.5 Mm and w = 0.3 Mm (top); and w for y0 = 0.5 Mm and A = 2 km s−1 (bottom).
In the text

	Fig. 11. Time–distance plots for Viy (left) and δTi/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 = y0 = 0.5 Mm.
In the text

	Fig. 12. Time–distance plots for δTi/T, collected along y for x = 0 Mm, for the pulse of its amplitude A = 2 km s−1, and width w = 0.3 Mm, launched from y = y0 = 0.5 Mm, in the case of By = 5 G (top left), By = 15 G (top right), By = 20 G (bottom left), and By = 25 G (bottom right).
In the text

	Fig. 13. Averaged over time and height relative to perturbed temperature of ions, H vs. By, for the pulse of its amplitude A = 2 km s−1 and width w = 0.3 Mm, launched from y = y0 = 0.5 Mm.
In the text