© ESO 2020

1. Introduction

The heating problem of the solar atmosphere is very challenging as it requires an explanation of the counter-intuitive 200−500 fold atmospheric temperature rise from about 5600 K at the solar surface over a height of about 2000 km. Two main families of heating mechanisms have been proposed: one uses microflaring while the other is based on energy carried by various waves present in the solar atmosphere. There is observational evidence for the existence of both microflares and waves in the solar atmosphere. The former likely contributes to the heating at higher altitudes, whereas the latter is expected to be responsible for the heating of lower atmospheric layers. The main unsolved heating problem refers to the mechanism by which the wave energy is generated and dissipated. A solution of this longstanding solar physics problem is proposed in this paper.

The most efficient source of waves is the convection zone. Lying below the photosphere, this zone is populated by highly turbulent motions that generate magnetohydrodynamic gravity waves (e.g., Musielak & Ulmschneider 2001). Recent observations reveal a diversity of waves in the solar atmosphere and they also demonstrate that some of these waves carry enough energy to heat this atmosphere (e.g., Jess et al. 2009; Srivastava et al. 2018). Propagation of these waves in the atmosphere is affected by cutoff frequencies. Acoustic-type waves with frequencies higher than their cutoffs are propagating, but waves with lower frequencies are evanescent and their amplitudes decay with height.

The concept of cutoff was originally introduced by Lamb (1909) and later modified to be applicable to the solar atmosphere by others (e.g., Roberts 2004; Musielak et al. 2006). Strong observational evidence for the existence of cutoffs in the atmosphere and their variations with the atmospheric height was presented by Wiśniewska et al. (2016) and Kayshap et al. (2018). Several theoretical attempts have been made to fit the available observational data (e.g., Murawski et al. 2016; Wójcik et al. 2018) without achieving satisfactory agreement. However, in the most recent study by Wójcik et al. (2019a), good agreement between the numerical simulations and observations was finally achieved, and the authors pointed out the importance of partially ionized effects that were accounted for by their numerical code. Similar code is used in this paper but with the exception that all nonideal and nonadiabatic terms are switched on. However, the main goal of this paper is different from that of Wójcik et al., as here we concentrate on wave generation and dissipation, and the solar atmosphere wave heating problem.

Waves in the solar atmosphere may form shocks that heat the plasma locally (e.g., Ulmschneider et al. 1978; Carlsson & Stein 1995). Fawzy et al. (2002) constructed theoretical models of the chromosphere that are based on the waves generated in the convection zone and using the dissipation of wave energy by shock formation (Hillier et al. 2016). The constructed models were purely theoretical, one-dimensional (1D) and they included both nonmagnetic and magnetic regions as well as full treatment of nonlocal thermodynamic equilibrium (NLTE) ionization and radiation. The only free parameter in these models was a filling factor, which was used to explain different levels of activity observed in solar-like stars. The models were used to compute profiles of Ca II and Mg II lines, which were compared to observational data. The comparison shows good agreement between the theoretically predicted line profiles and observations, and also demonstrates that the theoretically predicted temperature variations in the solar atmosphere are in an agreement with the semi-empirical VAL model (Vernazza et al. 1981). The main difference between the approach of Hillier et al., which was based on fully ionized plasma, and the one presented in this paper is that our numerical model is 2D and properly accounts for partially ionized plasma in the solar atmosphere.

In the chromosphere models described above, and in most of the previous theoretical studies, i.e., studies of Hillier et al. (2016), the effects of partially ionized plasma in the gravitationally stratified solar atmosphere were often neglected or considered under some assumptions. For instance, Wójcik et al. (2018) and Kuźma et al. (2019) devised 1D models of the dynamics of ions and neutrals and then the models were extended to 2D geometry (e.g., Srivastava et al. 2018; Wójcik et al. 2019a,b). Moreover, Maneva et al. (2017) developed a 2D model of nonadiabatic and nonideal atmosphere to show that ion magnetoacoustic-gravity and neutral acoustic-gravity waves (e.g., Vigeesh et al. 2017) lead to local heating of a magnetic flux tube. In addition, Martínez-Sykora et al. (2017) and Khomenko et al. (2018) included neutrals in their respectively 2D and 3D models of self-generated solar granulation. The contribution of neutrals was taken into account by amending a few extra terms in the induction equation. These models were recently generalized to dynamic ions and neutrals in the chromosphere by Popescu Braileanu et al. (2019) who reported the development of a new two-fluid code and its numerical verification. We present here an approach in which ions and neutrals are included self-consistently in our numerical model, and their dynamical effects on wave propagation and dissipation are fully accounted for. Our approach provides new insight into the wave heating problem.

Our paper is organized as follows: our two-fluid (ions + electrons and neutrals) model is described in Sect. 2; wave heating resulting from our numerical simulations is presented and discussed in Sect. 3; our conclusions are given in Sect. 4

2. Two-fluid model of a partially ionized solar atmosphere

We consider a gravitationally stratified and partially ionized solar atmosphere whose evolution is described by the set of two-fluid equations (e.g., Zaqarashvili et al. 2011; Leake et al. 2012; Martínez-Gómez et al. 2018; Soler et al. 2013; Oliver et al. 2016; Ballester et al. 2018, and references therein). We take into account hydrogen as a main dynamic plasma ingredient; the influence of heavier plasma species is considered with use of the OPAL solar abundance model (e.g., Vögler et al. 2004). Additionally, for simplicity we assume that ions and electrons consist of a single ion-electron fluid, while neutrals are described by a second fluid. Interaction between both fluids occurs through the ion–neutral collisions. In this model we neglect effects of recombination, ionization, and charge exchange, as well as effects of electrons exerted on ions, mimicked by extra terms in the generalized Ohm’s law (e.g., Khomenko 2017). As a result of complexity of the present model, these effects are left for future studies.

2.1. Two-fluid equations for ions and neutrals

We write the equations that describe dynamics of neutral and ion + electron fluids as

$$\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{i}}}{\partial t}+\mathrm{\nabla }·({\varrho }_{\mathrm{i}}{\mathit{V}}_{\mathrm{i}})=0,\hfill \end{array}$$(2)

$$\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}}\mathbf{I})\hfill \\ \hfill & ={\alpha }_c({\mathit{V}}_{\mathrm{i}}-{\mathit{V}}_{\mathrm{n}})+{\varrho }_{\mathrm{n}}\mathit{g},\hfill \end{array}$$(3)

$$\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}}\mathbf{I})\hfill \\ \hfill & =\frac{1}{\mu }(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}+{\alpha }_c({\mathit{V}}_{\mathrm{n}}-{\mathit{V}}_{\mathrm{i}})+{\varrho }_{\mathrm{i}}\mathit{g},\hfill \end{array}$$(4)

$$\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}$$(5)

$$\begin{array}{cc}& \frac{\partial E_{\mathrm{n}}}{\partial t}+\mathrm{\nabla }·[(E_{\mathrm{n}}+p_{\mathrm{n}}){\mathit{V}}_{\mathrm{n}}]\hfill \\ \hfill & ={\alpha }_c{\mathit{V}}_{\mathrm{n}}·({\mathit{V}}_{\mathrm{i}}-{\mathit{V}}_{\mathrm{n}})+Q_{\mathrm{n}}^{\mathrm{in}}+q_{\mathrm{n}}+{\varrho }_{\mathrm{n}}\mathit{g}·{\mathit{V}}_{\mathrm{n}},\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}·\mathit{B})]\hfill \\ \hfill & ={\alpha }_c{\mathit{V}}_{\mathrm{i}}·({\mathit{V}}_{\mathrm{n}}-{\mathit{V}}_{\mathrm{i}})+Q_{\mathrm{i}}^{\mathrm{in}}+Q_R^i+q_{\mathrm{i}}+{\varrho }_{\mathrm{i}}\mathit{g}·{\mathit{V}}_{\mathrm{i}},\hfill \end{array}$$(7)

where the heat production terms are

$$\begin{array}{c}\hfill Q_{\mathrm{n}}^{\mathrm{in}}={\alpha }_c(\frac{1}{2}{|\mathrm{\Delta }\mathit{V}|}^2+a\mathrm{\Delta }T),\end{array}$$(8)

$$\begin{array}{c}\hfill Q_{\mathrm{i}}^{\mathrm{in}}={\alpha }_c(\frac{1}{2}{|\mathrm{\Delta }\mathit{V}|}^2-a\mathrm{\Delta }T),\end{array}$$(9)

with ΔV = V_i − V_n, ΔT = T_i − T_n, and a = 3k_B/(m_H(μ_i+μ_n)). Here I is the identity matrix, and subscripts i, n, and e correspond respectively to ions, neutrals, and electrons. The symbols 𝜚_i, n denote mass densities, V_i, n velocities, p_ie, n ion + electron and neutral gas pressures, B magnetic field, and T_i, n are temperatures specified by ideal gas laws,

$$\begin{array}{c}\hfill p_{\mathrm{n}}=\frac{k_{\mathrm{B}}}{m_{\mathrm{H}}{\mu }_{\mathrm{n}}}{\varrho }_{\mathrm{n}}{T}_{\mathrm{n}},p_{\mathrm{ie}}=\frac{k_{\mathrm{B}}}{m_{\mathrm{H}}{\mu }_{\mathrm{i}}}{\varrho }_{\mathrm{i}}{T}_{\mathrm{i}},\end{array}$$(10)

where the gravity vector is g = [0, −g, 0] with its magnitude g = 274.78 m s⁻², α_c is the coefficient of collisions between ion and neutral particles, q_i, n describe thermal conduction (Spitzer 1962), $Q_{\mathrm{R}}^{\mathrm{i}}$ represents radiative losses, μ_ie = 0.58 and μ_n = 1.21 are the mean masses of electrons + ionized plasma species and neutrals, respectively, and they are taken from the OPAL solar abundance model (e.g., Vögler et al. 2004), m_H is the hydrogen mass, k_B is the Boltzmann constant, γ = 1.4 is the specific heats ratio, and μ is the magnetic permeability of the medium. All other symbols have their standard meaning.

The radiative losses term, $Q_{\mathrm{R}}^{\mathrm{i}}$, consists of two separate parts: thick cooling which operates in the low atmospheric layers and thin cooling, which works in the upper atmospheric regions. In every time-step we calculate the optical depth, starting from infinity and ending at the bottom boundary. For optical depths higher than 0.1 we use the thick cooling approximation of Abbett & Fisher (2012). For optical depths lower than 0.1 we use the thin cooling curve (Cox & Tucker 1969), which is interpolated by a tenth-order polynomial. In order to mimic energy transport into the simulation box from lower lying convective regions, we supplemented the radiative term by the heating term which overwhelms cooling by 5% at optical depths τ >  10.

Following Zaqarashvili et al. (2011), we assume coefficient of collisions between ions and neutrals α_in = α_ni = α_c, and to estimate its value we use the following formula provided by Braginskii (1965):

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

Here the collisional cross-section σ_in is taken from the quantum-mechanical model of Vranjes & Krstic (2013) who showed that the classical hard-sphere model may lead to underestimation of the cross-section values. For typical chromospheric plasma temperature in the range of 6 × 10³ − 10⁴ K this cross-section is equal to 1.89 × 10⁻¹⁸ m², that is about three orders of magnitude larger than in the hard-sphere model. We experienced for acoustic waves that the released heating due to ion-neutral collisions is higher for the classical cross-section than for its quantum analog (Kuźma et al. 2019). We note that the ion–neutral collision frequency differs from the neutral–ion collision frequency (Ballester et al. 2018), and they are given as ν_in = α_c/𝜚_i and ν_ni = α_c/𝜚_n, respectively.

2.2. Initial hydrostatic equilibrium of the solar atmosphere

Here, we consider the initial conditions of the solar atmosphere. We assume that the atmosphere is initially in hydrostatic (V_i, n = 0) equilibrium. The equilibrium mass density profiles of ions and neutrals (Fig. 1) are specified uniquely by the temperature profile, which is taken to be the same for ions and neutrals (e.g., Oliver et al. 2016) from the quiet solar atmosphere model of Avrett & Loeser (2008). We note that the ion mass density is much lower than the neutral mass density in the photosphere and the low chromosphere; they become comparable about 900 km below the transition region which is located at y ≈ 2.1 Mm. Higher up, the neutrals are less abundant than ions; in the corona the mass density of neutrals experiences a sudden fall-off with height. As a result of ion and neutral distributions, the ionization degree depends on height in the solar atmosphere; in the photosphere, where the temperature is only about 5800 K, the plasma is very weakly ionized with only about one ion for every thousand neutrals, while in the solar corona the plasma there is essentially fully ionized.

Fig. 1. Initial vertical profiles of ion, 𝜚i, (solid line) and neutral, 𝜚n, (dashed line) mass densities.

Interaction between ion + electron and neutral fluids occurs through the ion–neutral collisions that are determined by the ion–neutral friction coefficient (Eq. (11)). Vertical profile of α_c is illustrated in Fig. 2 (left). Below the photosphere, α_c retains its high value of 20 km s⁻¹, and therefore both fluids are strongly coupled there, acting as a single fluid as a result. At the bottom of the photosphere, that is at y = 0, the value of α_c begins to decrease with altitude. Higher up, both fluids decouple, enabling ion-neutral drift, ΔV = V_i − V_n, to take place in the chromospheric and transition region layers. This drift is essential for collisional heating (Eqs. (8) and (9)). The sound speed c_s can be given as

$$\begin{array}{c}\hfill c_{\mathrm{s}}(y)=\sqrt{\frac{\gamma (p_{\mathrm{h}\mathrm{ie}}(y)+p_{\mathrm{h}\mathrm{n}}(y))}{{\varrho }_{\mathrm{h}\mathrm{i}}(y)+{\varrho }_{\mathrm{h}\mathrm{n}}(y)}}·\end{array}$$(12)

Fig. 2. Vertical profiles of ion-neutral friction coefficient, αc, (left panel) and sound speed, cs, (right panel) in units of g cm−3 s−1 and km s−1.

In the photosphere, c_s is close to 6 km s⁻¹ (Fig. 2, right). Immediately below the transition region, c_s(y = 2 Mm) ≈ 11 km s⁻¹. Higher up, c_s grows, following the trend of $\sqrt{T(y)}$; this growth is abrupt immediately above the transition region, and at larger altitudes c_s rises slowly to a value of about 150 km s⁻¹ at y ≈ 30 Mm (not shown).

We adopt here the uniform magnetic field with its vertical and transverse components given as B = [0, B₀, B₀] with $B_0=15/\sqrt{2}$ Gs. A nonzero transversal magnetic field permits transversal Alfvén waves to be present in the system. These waves are linearly coupled to magnetoacoustic waves and are essential in the model; as in the nonlinear regime they are capable of driving vertical flow (e.g., Hollweg 1986; Murawski 1992; Shestov et al. 2017). The mode conversion between slow-acoustic and Alfvén modes in regions of β ≈ 1 can take place in this system. A nonlinear coupling between slow and Alfvén waves near these regions has been studied by Zaqarashvili & Roberts (2006). A detailed discussion on two-fluid Alfvén waves can be found in Kuźma et al. (in prep.). The vertical profile of the corresponding Alfvén speed, $c_{\mathrm{A}}(x,y)=B/\sqrt{\mu ({\varrho }_{\mathrm{i}}+{\varrho }_{\mathrm{n}})}$, is illustrated in Fig. 3. In the upper photosphere and lower chromosphere, c_A attains values lower than 10 km s⁻¹ up to y ≈ 1 Mm. Higher-up, as a result of low ion mass density 𝜚_i, c_A rises until the transition region, where it experiences a sudden growth to values close to 800 km s⁻¹.

Fig. 3. Profile of the initial Alfvén speed, cA, versus height y.

3. Numerical simulations of two-fluid wave heating

We use JOANNA code (Wójcik et al. 2018) to solve the two-fluid equations numerically. We set in our numerical experiments the Courant-Friedrichs-Lewy (Courant et al. 1928) number equal to 0.3 and choose a second-order accuracy in space and a second-order accurate Runge-Kutta method (Durran 2010) for integration in time, supplemented by adopting the Harten-Lax-van Leer discontinuities (HLLD) approximate Riemann solver (Miyoshi & Kusano 2005) and the divergence of magnetic field cleaning method of Dedner et al. (2002).

The simulation box is specified along horizontal (x-) and vertical (y-) directions as ( − 5.12 <  x <  5.12) Mm × ( − 5.12 <  y <  30) Mm. This allows us to simulate a convectively unstable region just below the photosphere which occupies the layer 0 ≤ y ≤ 0.5 Mm. Below the height y = 5.12 Mm, we set a uniform grid with cell size 20 km × 20 km, while higher up at the base of the corona, which is located at y = 2.1 Mm, we stretch the grid along the y-direction dividing it into cells whose size grows with height. Two different types of boundary conditions are implemented at four edges of the simulation box. All plasma quantities at the top and bottom boundaries are set to their hydrostatic values. At the right and left boundary we use periodic boundary conditions. With the use of the stretched grid above y = 5.12 Mm we damp any incoming signal in layers close to the upper boundary. We find that such a stretched grid, supplemented by fixing all plasma quantities to their hydrostatic conditions at the top boundary, significantly reduces spurious reflections of the incoming signal. Additionally, we implemented a plasma inflow at the bottom boundary with its vertical velocity equal to 0.3 km s⁻¹. This allows us to compensate the outflowing plasma mass losses. The value of inflow velocity is calculated by assuming that the incoming flux compensates the radiative losses (Chatterjee 2020).

In our numerical simulations we consider plasma that is initially (at t = 0 s) set to its hydrostatic state with the vertical and transverse components of magnetic field equal to $10/\sqrt{2}$ Gs. Later on, horizontal motions of plasma due to spontaneously generated granulation reshuffle magnetic field lines in the photosphere in the form of small flux tubes (Fig. 4, solid lines). These horizontal plasma motions reveal mean velocity close to 2.5 km s⁻¹ at y = 0 Mm (not shown) with a maximum value of up to 5 km s⁻¹, and thus they remain subsonic. This figure also reveals the presence of a well-defined cool and dynamic chromosphere above the photosphere. This chromosphere layer is filled with both hot and cold jets visible at x ≈ 2.5 Mm.

Fig. 4. Spatial profile of logarithm of ion temperature, Ti, expressed in units of Kelvins (colormap) overlayed by magnetic field lines at t = 5000 s.

Hot plasma associated with granulation (Fig. 4, colormap) drifts towards the chromosphere and forms plasma upflows (Fig. 5, top). At the bottom of the chromosphere this plasma radiates its energy and gets cooled. Between sibling granules the downdrafts with cold plasma flowing down settle in the convection zone (Fig. 5, top). Both these upflows and downdrafts are accompanied by two-fluid waves (Fig. 5, bottom), which in our model atmosphere are neutral acoustic-gravity, ion magnetoacoustic-gravity (Vigeesh et al. 2017), and Alfvén waves.

Fig. 5. Spatial profile of the vertical component of ion velocity, Vi, expressed in units of 1 km s−1 (top) and its spatially averaged Fourier power periods vs. height (bottom) at t = 5000 s.

The periods of the excited waves, which we obtained from the Fourier spectra for ion vertical velocity, are plotted in Fig. 5 (bottom). We note that at the height within the range of 0 Mm <  y <  0.5 Mm, the dominant wave period P is about 340 s. However, higher up, for y >  0.5 Mm, the signal associated with large values of P decays fast with y, and the main wave period is associated with P ≈ 180 s. It was demonstrated that acoustic waves of P ≈ 300 s are evanescent in the photosphere (e.g., Wójcik et al. 2018) but they are filtered out in higher layers. Thus, the package of waves consists of the main wave-period of about 180 s (Fig. 5, bottom). The described filtering process dominates in the chromosphere and it lowers periods and allows the waves to reach the transition region located at y = 2.1 Mm, where the waves become partially reflected back to the chromosphere. As a result of ion–neutral collisions, the energy of these excited waves is dissipated. This dissipation is most effective for largest dispatches between ion and neutral velocities, and these waves may convert their energy into heat mostly in the chromosphere (Fig. 6), compensating radiative and thermal losses.

Fig. 6. Spatial profile of the absolute difference between ion and neutral velocities (|Vi − Vn|) expressed in units of 1 km s−1 at t = 5000 s.

A strong damping of low-frequency (with period > 80 s) Alfvén waves does not follow the conclusions drawn from the analytical studies of Zaqarashvili et al. (2013) who showed that waves with short wave-periods are rapidly damped by ion–neutral collisions. However, our two-fluid model describes highly dynamic solar atmosphere (that nevertheless remains in quasi-equilibrium) and differs significantly from those developed so far. A signal present in V_iz is partially converted through ion-neutral collisions into a signal in V_nz. It follows from the dispersion relation that V_nz is associated with a purely imaginary frequency which corresponds to a strongly damped wave (Zaqarashvili et al. 2013). This effect along with the energy transfer from V_iy into V_nz causes Alfvén waves, which correspond to V_iz, to be strongly damped. We found that longer wave periods are not the main cause of the heating of the solar atmosphere, which agrees with the numerical results of Kuźma et al. (2019) who found that for acoustic waves, with wave periods longer than 80 s, the plasma heating decreases.

Figure 7 shows the vertical profile of horizontally averaged ionization, 𝜚_i/𝜚_n, at three moments of time, namely t = 0 s (solid line), that is, initial ionization, t = 2000 s (dashed line), and t = 4000 s (dotted line), which corresponds to fully developed and ongoing granulation. From this plot we can see, that 𝜚_i/𝜚_n drops by about one order of magnitude in the upper chromosphere. The ionization levels in the photosphere, lower chromosphere, and corona remain close to their initial values.

Fig. 7. Horizontally averaged vertical profile of ionization fraction, 𝜚i/𝜚n, at t = 0 s (solid line), t = 2000 s (dashed line) and t = 4000 s (dash-dotted line).

The numerically obtained quasi-stationary solar atmosphere is the result of the subtle balance between radiative cooling, thermal conduction, and collisional energy dissipation mechanism. This equilibrium is illustrated in Fig. 8 in the form of horizontally averaged collisional heating and negative cooling rates versus height. We note that by definition, the heating rate attains positive values, while cooling negative values. Deep below the photosphere, heating from the deeper layers overcomes cooling, and leads to hot plasma upflows in the convectively unstable layer. At the photosphere, that is for 0 <  y <  0.5 Mm, the optical opacity abruptly falls off and hot plasma rapidly radiates its energy. Thereafter, the heating term overwhelms the cooling term in this layer. Higher up, that is, in the chromosphere and low corona, these two terms reach a similar value, with the cooling term being slightly dominant in the corona. It is demonstrated that the wave energy dissipation, which is supplemented by viscous and ohmic heatings, is sufficient to balance radiative and thermal energy losses above the photosphere and to sustain a quasi-stationary solar atmosphere.

Fig. 8. Horizontally averaged vertical profiles of collisional heating (dotted line) and minus cooling term (solid line) vs. height y, drawn at t = 5000 s.

Figure 9 shows the temporal evolution of horizontally averaged T_i. We see that, as a result of cooling–heating balance in the chromosphere, the system arrives at a quasi-equilibrium state of the solar atmosphere. The chromosphere in our simulations is a highly dynamic medium filled with jets and spicules that transport mass to balance losses due to plasma outflows (Wójcik et al. 2019b). The well-defined transition region oscillates close to its initial position at y = 2.1 Mm.

Fig. 9. Time distance plot of the horizontally averaged Ti.

Figure 10 reveals that the radiative cooling and thermal energy losses are in equilibrium with collisional heating, sustaining a quasi-stationary atmosphere with the transition region located around y = 2.1 Mm. As a result of this balance, the numerically obtained vertical profile of ion temperature is close to the semi-empirical data of Avrett & Loeser (2008). A better fit would require tuning the inflow or the extra heating which we impose below the photosphere. However, such parametric studies are very costly as they require a number of jobs to be run. It must be pointed out that this agreement was obtained without invoking shocks, which were previously considered as the main source of local heating in the solar chromosphere (e.g., Ulmschneider et al. 1978; Carlsson & Stein 1995; Fawzy et al. 2002). Such shocks are occasionally generated in the upper chromosphere. However, our results show that the efficient wave heating by ion-neutral collisions prevents shock formation in the lower and middle chromosphere. Thus, a new insight into the wave heating problem is our demonstration that the waves dissipate their energy without forming shocks; another important aspect of this result is that short-period waves, which were the main source of heating of the lower chromosphere in the latter mentioned study, can now reach much higher layers of the chromosphere, including the transition region and corona. Our results allow us to obtain an average temperature distribution that is very similar to that in the semi-empirical model of Avrett & Loeser (2008), which is good confirmation of the validity of our approach and the code used to perform the simulations. Moreover, we also demonstrate that the ion temperature varies locally in the chromosphere, and the predicted range is consistent with the results of the numerical simulations reported by Carlsson & Stein (1995). The latter study shows that temperature variation grows with height, attaining the greatest values in the low corona. Figure 10 presents a similar trend, although the temperature variations are smaller than those of Carlsson & Stein (1995), which may result from the presence of ions and neutrals in the model we developed.

Fig. 10. Numerically obtained minimum Ti(y) (blue dashed line) and maximum Ti(y) (red dashed line) at x = 0, and semi-empirical data of T(y) of Avrett & Loeser (2008; dotted line).

Recent observations of Grant et al. (2018) showed that Alfvén waves are steepening into shocks at lower heights of the solar atmosphere. With observed local temperature enhancements of 5% it is the first observational evidence of Alfvén waves heating chromospheric plasma in a sunspot umbra through the formation of shock fronts. In our simulations, thermal energy is released during collisions between ions and neutrals which may lead to ion–neutral velocity drift, and this is a pure two-fluid effect. In the case of a shock, kinetic energy at the front of this shock is converted to thermal energy as a result of the friction between ion + electron fluid and neutral species. This phenomenon takes place not only at a shock but also at all other waves in general. However, it is more efficient for steeper wave profiles and its impact on overlaying plasma temperature depends on ion–neutral drift. In conclusion, it is therefore largest in the case of shocks.

Finally, it must also be emphasized that similar code was used by Wójcik et al. (2019a,b) to propose that the solar wind origin is associated with granulation, and to theoretically predict the variations of acoustic cutoff frequency in the solar atmosphere, which are in good agreement with the observational results obtained by Wiśniewska et al. (2016) and Kayshap et al. (2018). This is another independent verification of the validity of our approach and the code, and it shows that the effects of partial ionization become very important in lower and middle chromospheric layers.

4. Conclusions

We performed numerical simulations of two-fluid waves in a partially ionized solar atmosphere with radiation and thermal conduction taken self-consistently into account. The considered neutral acoustic-gravity, Alfvén, and ion magnetoacoustic gravity waves were generated by spontaneously evolving and self-organizing convection with granulation cells. We found that the granulation excites a wide spectrum of wave periods (Fig. 5, bottom) and that only sufficiently short-period acoustic and long-period gravity waves propagate through the atmosphere and reach the corona (Vigeesh et al. 2017).

Part of the energy carried by these waves is dissipated in the photosphere and chromosphere because of ion–neutral collisions. As a result, these two-fluid waves generated by the granulation effectively heat the plasma, compensating radiative and thermal energy losses there. The atmosphere also gets extra heat from viscosity and magnetic diffusivity, which are most effective at the places of sheared flow and oppositely orientated magnetic field lines. This heating takes place at the downdrafts below the photosphere, contributing to plasma temperatures there and leading to a quasi-stationary solar atmosphere with its temperature distribution being close to the semi-empirical model of Avrett & Loeser (2008).

This is the first time that agreement with the observationally based semi-empirical model has been obtained without involving heating by shocks. We conclude that our results unravel the main mechanism of the wave heating quiet regions of the chromosphere, and that the basic physical process of this mechanism is collisions between neutrals and ions. Our results have important implications for the construction of theoretical solar atmospheric models in the near future as they show that the effects of partially ionized plasma must be included in these models and that they may prevent shock formation, which was a dominant physical process for waves in a fully ionized solar atmosphere.

Acknowledgments

The JOANNA code was developed by Darek Wójcik at UMCS. This work was done within the framework of the projects from the Polish Science Center (NCN) Grant Nos. 2017/25/B/ST9/00506 and 2017/27/N/ST9/01798.

References

All Figures

	Fig. 1. Initial vertical profiles of ion, 𝜚i, (solid line) and neutral, 𝜚n, (dashed line) mass densities.
In the text

	Fig. 2. Vertical profiles of ion-neutral friction coefficient, αc, (left panel) and sound speed, cs, (right panel) in units of g cm−3 s−1 and km s−1.
In the text

	Fig. 3. Profile of the initial Alfvén speed, cA, versus height y.
In the text

	Fig. 4. Spatial profile of logarithm of ion temperature, Ti, expressed in units of Kelvins (colormap) overlayed by magnetic field lines at t = 5000 s.
In the text

	Fig. 5. Spatial profile of the vertical component of ion velocity, Vi, expressed in units of 1 km s−1 (top) and its spatially averaged Fourier power periods vs. height (bottom) at t = 5000 s.
In the text

	Fig. 6. Spatial profile of the absolute difference between ion and neutral velocities (|Vi − Vn|) expressed in units of 1 km s−1 at t = 5000 s.
In the text

	Fig. 7. Horizontally averaged vertical profile of ionization fraction, 𝜚i/𝜚n, at t = 0 s (solid line), t = 2000 s (dashed line) and t = 4000 s (dash-dotted line).
In the text

	Fig. 8. Horizontally averaged vertical profiles of collisional heating (dotted line) and minus cooling term (solid line) vs. height y, drawn at t = 5000 s.
In the text

	Fig. 9. Time distance plot of the horizontally averaged Ti.
In the text

	Fig. 10. Numerically obtained minimum Ti(y) (blue dashed line) and maximum Ti(y) (red dashed line) at x = 0, and semi-empirical data of T(y) of Avrett & Loeser (2008; dotted line).
In the text