© ESO 2021

1. Introduction

One of the main problems of heliophysics concerns the source of the thermal energy required to heat the solar atmosphere, the lowest layer of which is about 500 km thick and is called the photosphere. The chromosphere, which is about 1.5 Mm thick, caps this layer. A narrow transition region separates the chromosphere from the solar corona above it. Temperatures in the photosphere vary from about 5600 K at its bottom to about 4300 K at its top. The temperature of the middle of the chromosphere rises to about 10⁴ K (±800 − ±1800 km above the solar surface). Because of the temperatures in the photosphere and low and middle chromosphere, the plasma in these layers is only partially ionized. Recombination, photoionization, and chemical ionization were discussed by Carlsson & Stein (2002). At the top of the photosphere, the ionization degree is as low as 10⁻⁴, and so there are about 10⁴ neutrals per ion. In the chromosphere, the ratio of ions to neutrals grows with height, and therefore with temperature (Khomenko et al. 2014). As a result of its 1 − 3 MK average temperature, the plasma in the corona is almost fully ionized.

The plasma density is also higher in the chromosphere and thus it radiates more than the corona. Hence, an additional mechanism of heating is required to compensate for these radiative losses (Narain & Ulmschneider 1996). One of the sources of this heating may be associated with collisions between particles, which lead to wave damping (Erdélyi & James 2004). The idea of heating the solar chromosphere by acoustic waves is very old and was first developed by Biermann (1946) and Schwarzschild (1948), who suggested that acoustic waves play a major role in the heating of this atmospheric layer. Later, the problem was widely studied, for example by Carlsson & Stein (1995), Ulmschneider et al. (2003), and Kuźma et al. (2019). The latter authors focused on the acoustic waves with wave periods within the range 30 − 200 s and showed that these waves are able to heat the chromosphere. The idea of heating the solar atmosphere without shock waves was studied by Wójcik et al. (2020) who showed that the thermal energy release ratio during ion–neutral collisions is sufficient to balance the radiative and thermal energy losses. Recently, Kuźma et al. (2021) focused on the studies of ion magnetoacoustic and neutral acoustic waves in a magnetic arcade and found that, in the chromosphere, the wave characteristics strongly depend on the magnetic field configuration. The advantages of using the two-fluid model over the single-fluid model were studied by Zaqarashvili et al. (2011b). The single-fluid (MHD) model is sufficiently accurate for slow processes but fails on timescales shorter than the ion–neutral collision time. Hence, the two-fluid approximation should be used for studying faster processes.

The problem of atmospheric heating was also studied in the context of Alfvén waves (Piddington 1956; Osterbrock 1961). The mechanisms by which these waves are dampened by ion–neutral collision were studied by De Pontieu et al. (2001). In the lower chromosphere, collisional effects were found to be the most important component of the dissipation mechanism (Khodachenko et al. 2004). As a matter of fact, Alfvén waves can dissipate their energy in the chromosphere as a result of ion–neutral collisions (Leake et al. 2005; Goodman 2011; Tu & Song 2013; Zaqarashvili et al. 2013; Arber et al. 2016; Shelyag et al. 2016; Soler et al. 2017). However, this mechanism is only effective in the chromosphere and not in the photosphere or in the corona because of weak magnetization and strong ionization, respectively. In other studies, Zaqarashvili et al. (2013) and Soler et al. (2017) focused on the dissipation of torsional Alfvén waves and showed that they can by effectively dissipated in the chromosphere.

The mechanisms considered in the present paper for solar coronal heating are essentially based on magnetic field reconnection and resonant dissipation and/or phase-mixing of magnetohydrodynamic (MHD) waves. The latter mechanism was reviewed recently by Van Doorsselaere et al. (2020), who showed that various wave modes may contribute to the heating of the solar corona. The effect of the cooling and heating on the dynamics of the MHD waves was investigated by Nakariakov et al. (2017) who concluded that the perturbations of the thermal equilibrium, for example by slow magnetoacoustic waves, lead to a heating–cooling unbalance which acts as an energy exchange mechanism. For more recent studies of thermal unbalance, see also Prasad et al. (2021) who studied the damping of slow magnetosonic waves based on a new dispersion relation. The problem of damping slow magnetosonic waves in the solar corona was also investigated by Duckenfield et al. (2021), considering nonzero β plasma, and Kolotkov et al. (2019) who presented three potential regimes of wave evolution depending on characteristic timescales of the thermal unbalance. Recently, Kolotkov et al. (2020) constrained the coronal heating function by observations of slow magnetoacoustic waves.

The solar corona extends smoothly into the solar wind. Studies of the solar wind were initiated by Biermann (1951) and Parker (1965). Parker (1958) postulated that the interplanetary space is not void but that it is filled by the solar wind. An early model of the solar wind with the power spectrum of Alfvén and kink waves was developed by Tu (1987). Later on, a lot of attention was paid to the study of the formation of the solar wind by plasma outflows. For instance, Tu et al. (2005) found these outflows in coronal funnels at altitudes between 5 and 20 Mm above the photosphere. Later, the problem was reconsidered by Wójcik et al. (2020) who proposed that jets and associated plasma outflows generated by the solar granulation may result in the formation of the solar wind higher up. Moving upwardly, the plasma requires a source of momentum which can be provided by MHD waves (Ofman 2005; Marsch 2006; De Pontieu et al. 2007). Alfvén waves, for instance, may drive plasma outflows (Hollweg et al. 1982; Kudoh & Shibata 1999; Matsumoto & Suzuki 2012). A model for coronal heating and solar wind acceleration by low-frequency Alfvén wave turbulence was presented by van der Holst et al. (2014).

The goal of the present paper is to focus on magnetoacoustic waves in the context of heating the chromosphere and the associated plasma outflows, and therefore we leave the role of Alfvén waves for future studies. The novelty of this work is that we consider impulsively generated two-fluid magnetoacoustic waves in the context of solar chromospheric heating and plasma outflows. This paper is organized as follows. In Sect. 2, we present the two-fluid equations that are considered for describing the dynamics of the partially ionized solar atmosphere. Numerical results are presented in Sect. 3. In Sect. 4 we conclude this paper by summarizing the main results and briefly discussing them.

2. Two-fluid numerical model of the partially ionized solar atmosphere

2.1. Two-fluid equations

To describe the low layers of the solar atmosphere we use two-fluid equations for ions + electrons (_i, e) and neutrals (_n) treated as separate interacting fluids (Zaqarashvili et al. 2011a; Leake et al. 2014; Oliver et al. 2016; Maneva et al. 2017; Popescu Braileanu et al. 2019):

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

$$\begin{array}{cc}\hfill & \frac{\partial {\varrho }_n}{\partial t}+\mathrm{\nabla }·({\varrho }_n{\mathbf{V}}_n)=0,\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill & \frac{\partial ({\varrho }_i{\mathbf{V}}_i)}{\partial t}+\mathrm{\nabla }·({\varrho }_i{\mathbf{V}}_i{\mathbf{V}}_i+p_{\mathit{ie}}\mathbf{I})=\hfill \\ \hfill & {\varrho }_i\mathbf{g}+\frac{1}{\mu }(\mathrm{\nabla }\times \mathbf{B})\times \mathbf{B}+{\mathbf{S}}_m,\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill & \frac{\partial ({\varrho }_n{\mathbf{V}}_n)}{\partial t}+\mathrm{\nabla }·({\varrho }_n{\mathbf{V}}_{\mathrm{n}}{\mathbf{V}}_n+p_n\mathbf{I})={\varrho }_n\mathbf{g}-{\mathbf{S}}_m,\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill & \frac{\partial E_i}{\partial t}+\mathrm{\nabla }·[(E_i+p_{\mathit{ie}}+\frac{{\mathbf{B}}^{\mathbf{2}}}{2\mu }){\mathbf{V}}_i-\frac{\mathbf{B}}{\mu }({\mathbf{V}}_i·\mathbf{B})]=\hfill \\ \hfill & ({\varrho }_i\mathbf{g}+{\mathbf{S}}_m)·{\mathbf{V}}_i+Q_i,\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill & \frac{\partial E_n}{\partial t}+\mathrm{\nabla }·[(E_n+p_n){\mathbf{V}}_n]=({\varrho }_n\mathbf{g}-{\mathbf{S}}_m)·{\mathbf{V}}_n+Q_n,\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill & \frac{\partial \mathbf{B}}{\partial t}=\mathrm{\nabla }\times ({\mathbf{V}}_i\times \mathbf{B}),\mathrm{\nabla }·\mathbf{B}=0,\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill & E_i=\frac{{\varrho }_i{\mathbf{V}}_i^2}{2}+\frac{p_{\mathit{ie}}}{\gamma -1}+\frac{{\mathbf{B}}^2}{2\mu },E_n=\frac{{\varrho }_n{\mathbf{V}}_n^2}{2}+\frac{p_n}{\gamma -1}·\hfill \end{array}$$(8)

Here, I is the identity matrix, p_ie, n denote the ion + electron and neutral gas pressures, while 𝜚_i, n denote ion and neutral mass densities, and V_i, n the ion and neutral velocities. Moreover, g = [0, −g, 0] is gravitational acceleration with its magnitude g = 274.78 m s⁻², B is the magnetic field pointing in the y−direction, and μ the magnetic permeability of the medium. In Eqs. (3)–(6), S_m and Q_i, n correspond to the collisional momentum and energy exchange terms, respectively. They are given as (Meier & Shumlak 2012; Oliver et al. 2016)

$$\begin{array}{cc}& {\mathbf{S}}_m={\varrho }_i{\nu }_{\mathit{in}}({\mathbf{V}}_n-{\mathbf{V}}_i),\hfill \end{array}$$(9)

$$\begin{array}{cc}& Q_i=\frac{1}{2}{\nu }_{\mathit{in}}{\varrho }_i{({\mathbf{V}}_i-{\mathbf{V}}_n)}^2+\frac{k_B{\nu }_{\mathit{in}}{\varrho }_i}{(\gamma -1)m_n}(T_n-T_i),\hfill \end{array}$$(10)

$$\begin{array}{cc}& Q_n=\frac{1}{2}{\nu }_{\mathit{in}}{\varrho }_i{({\mathbf{V}}_i-{\mathbf{V}}_n)}^2+\frac{k_B{\nu }_{\mathit{in}}{\varrho }_i}{(\gamma -1)m_n}(T_i-T_n).\hfill \end{array}$$(11)

The ion–neutral collision frequency, ν_in, is defined as (Braginskii 1965; Ballester et al. 2018)

$$\begin{array}{c}\hfill {\nu }_{\mathit{in}}=\frac{4}{3}\frac{{\sigma }_{\mathit{in}}{\varrho }_n}{m_H({\mu }_i+{\mu }_n)}\sqrt{\frac{8k_B}{\pi m_H}(\frac{T_i}{{\mu }_i}+\frac{T_n}{{\mu }_n})}.\end{array}$$(12)

In the above equation, m_H is the hydrogen mass, σ_in = 1.4 · 10⁻¹⁵ cm² is the collision cross-section (Vranjes & Krstic 2013), μ_i = 0.58 and μ_n = 1.21 are the mean masses of ions and neutrals, respectively, and T_i, n are the temperatures specified by the ideal gas laws,

$$\begin{array}{c}\hfill p_{\mathit{ie}}=\frac{k_B}{m_i}{\varrho }_i{T}_i,\text{and}{p}_n=\frac{k_B}{m_n}{\varrho }_n{T}_n.\end{array}$$(13)

Here, m_i  =  m_Hμ_i, and m_n  =  m_Hμ_n, where m_H corresponds to the mass of a hydrogen atom, k_B is the Boltzmann constant, and γ  =  1.4 corresponds to the specific heat ratio. As the photospheric plasma contains bi-atomic molecules, we decided to consider γ = 1.4 instead of 5/3. The latter value is suitable for plasma that only consists of atoms. For simplicity, we neglected thermal conduction and other nonadiabatic and nonideal terms in the two-fluid equations. We should expect that thermal conduction will impact wave propagation by diffusing thermal energy into less localized and thus less intensively plasma-heated regions.

The two-fluid equations do not describe the top of the convection zone and the photosphere very well, which are included in the simulation domain, but they do accurately describe the chromosphere and the low corona, which are the focus of our simulation model. Problems with the two-fluid equations result from comparable values of characteristic gyro- and collision frequencies in the photosphere Khomenko et al. (2014), Fig. 1, top), while fluid models require that the collision frequencies be much higher than the gyro frequencies.

2.2. Magnetohydrostatic equilibrium

As an initial condition, we assume that the atmosphere is static (V_i, n = 0) and is in hydrostatic equilibrium with the ion and neutral pressure gradients being balanced by the gravity forces on both species:

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

Using the ideal gas laws in Eq. (13) and the y-components of Eq. (14), we express the background gas pressures and mass densities as

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

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

where Λ_i, n = k_BT/(μ_i, ng) are ion and neutral pressure scale heights, respectively, while p_i0 = 0.01 Pa and p_n0 = 0.0003 Pa denote the corresponding gas pressures at the reference level, taken here at y = y_r = 50 Mm.

We consider the height-dependent temperature, T(y), model of Avrett & Loeser (2008); see Fig. 1 (top). This model also determines the equilibrium mass densities and gas pressure profiles. As the electron mass is very small in comparison to the mass of neutrals and ions, we neglect the electron dynamics. We note that the ion mass density is lower than that of the neutral component in the photosphere and low chromosphere (Fig. 1, middle). However, above the altitude of y ≈ 0.9 Mm and especially above the transition region, which is located at y ≈ 2.1 Mm, neutrals are much less numerous than ions. The presence of gravity in the system leads to the appearance of an acoustic cut-off period (Fig. 1, bottom). The formula is given as

$$\begin{array}{c}\hfill P_{\mathrm{ac}}=\frac{4\pi \mathrm{\Lambda }}{C_S\sqrt{1+2\frac{\mathrm{d}\mathrm{\Lambda }}{\mathrm{d}y}}},\end{array}$$(17)

Fig. 1. Variation with height of the equilibrium temperature (top), mass densities (middle) of neutrals (dashed line) and ions (dotted line), and acoustic cut-off period (bottom).

where Λ is the pressure scale height, and $C_s=\sqrt{\gamma (p_i+p_n)/({\varrho }_i+{\varrho }_n)}$ is the sound speed.

The pressure scale height Λ is responsible for the propagation of the waves when their period P is smaller than the P_ac. Otherwise, if this value is greater than or equal to P_ac, the waves are evanescent (Horace 1909). It is worth noting that the P_ac considered here is a local quantity that varies with height (Musielak et al. 2006).

The hydrostatic equilibrium described above is overlaid by a uniform vertical magnetic field of magnitude B_y = 10 Gs. Such a magnetic field is current-free (∇ × B/μ = 0) and therefore force-free ((∇ × B)×B/μ = 0), and thus does not alter the hydrostatic equilibrium.

2.3. Impulsive perturbations

At t = 0 s, when the numerical simulations are started, we perturb the magnetohydrostatic equilibrium with a localized signal in the vertical components of the ion and neutral velocities:

$$\begin{array}{c}\hfill V_{\mathit{iy}}(y,t=0)=V_{\mathit{ny}}(y,t=0)=Aexp(-{\left(\frac{y-y_0}{w}\right)}^2),\end{array}$$(18)

where A is the amplitude of the pulse, w its width, and y₀ the height from which the pulse is launched. We consider the following cases: A = 0.5 km s⁻¹, A = 1 km s⁻¹, A = 2 km s⁻¹, and A = −5 km s⁻¹, combined with different widths of the pulse, namely w = 0.1 Mm, w = 0.25 Mm, and w = 0.3 Mm; and various launching heights, namely y₀ = 0 Mm, y₀ = 0.25 Mm, and y₀ = 0.5 Mm. The first two values of the amplitude A represent a typical flow associated with the solar granulation (e.g., Falco et al. 2017), while the third value corresponds to a stronger signal, and the last value mimics downdrafts. The implied value of the pulse amplitude A = 0.5 km s⁻¹ is about six times larger than the local Alfvén speed, which for the value of B_y = 10 Gs used here is about C_A = 0.08 km s⁻¹ at y = 0 Mm and much smaller than the local sound speed which is about C_S = 7.5 km s⁻¹ at y = 0 Mm. The characteristic width of the pulse, w = 300 km, is comparable to the size of the smallest granules. The difference between the positive and negative values of the amplitude results in the ion and neutral fluids being pushed upwards and downwards, respectively. However, both negative and positive pulses spread into upward and downward propagating waves. We note that y₀ = 0 Mm corresponds to the bottom of the photosphere, y₀ = 0.25 Mm to the middle of the photosphere, and y₀ = 0.5 Mm is associated with its top. For simplicity, we consider only pulses in the ion and neutral vertical velocity components. The physical nature of realistic perturbations is indeed more complex than the dependence presented by Eq. (18). Moreover, we consider an impulsive perturbation launched at t = 0 s, without first relaxing the system, because the relative numerical errors are small in all atmospheric layers (Fig. 2).

Fig. 2. Time–distance plots for Viy/Cs (top) and δTi/T (bottom) in the case of no pulse (A = 0).

3. Numerical results

The one-dimensional numerical simulations are performed with the JOANNA code (Wójcik et al. 2020), which solves the above-mentioned two-fluid equations. We assume that the system has a slab geometry and that it is invariant along the x- and z-directions (∂/∂x = ∂/∂z = 0), and the z−components of the velocities and the magnetic field are both set to zero (V_iz = V_nz = B_z = 0). Applying these simplifications, we remove Alfvén waves from the system which guides (affected by gravity) magnetoacoustic waves. Internal-gravity waves are absent in the system as they are unable to propagate along the gravity action.

The simulation domain is specified as −0.5 ≤ y ≤ 60.0 Mm. The distance below 6.87 Mm is divided into 1024 cells, and therefore the size of the grid cells is Δy ≈ 7.2 km. Higher up, the grid is stretched and the remainder of the domain is covered by 128 cells. At the top and bottom boundaries of the simulation region, we set and hold all plasma quantities fixed to their equilibrium values, given by Eqs. (15) and (16).

3.1. Numerical test

Figure 2 illustrates time–distance plots for the ratio between the vertical component of the ion velocity and the local sound speed, V_iy/C_s, (top) and the relative perturbed temperature of ions, δT_i/T = (T_i − T)/T, (bottom) in the absence of any initial pulse (i.e., for A = 0 in Eq. (18)). The signals present in both panels result from numerical errors which are inherent features of every numerical code. As the max |V_iy/C_s|≈0.08 occurs in the corona, and it indicates that the vertical velocities resulting from the discretization are small in comparison to the local sound speed, we infer that the numerical errors induced in the velocity component V_iy are negligible. In the lower regions of the atmosphere, the numerical errors decrease to max |V_iy/C_s|≈0.02 in the chromosphere, and max |V_iy/C_s|≈0.0002 in the photosphere. Similarly, max |δT_i/T|≈0.065, showing that the numerical errors detected in the perturbed ion temperature are also negligibly small. In contrast to the numerical errors in velocity, the errors in the perturbed temperature decrease with height to max |δT_i/T|≈0.03 in the chromosphere, and max |δT_i/T|≈0.02 in the corona. The apparent periodicity in δT_i/T (Fig. 2, bottom) may coincide with the acoustic cut-off period of the numerical noise.

3.2. Magnetoacoustic waves

3.2.1. Pulse amplitude effects

We first considered an initial pulse at the bottom of the photosphere (y₀ = 0 Mm). The initial pulse, given by Eq. (18), excites magnetoacoustic waves. These waves can get shocked in the chromosphere (Snow & Hillier 2021). The formation of shock waves is clearly seen in the time–distance plots for V_iy (Fig. 3, top), and it is stronger for a larger value of |A| (right-top). These shocks heat plasma through ion–neutral collisions. Similarly to the gravity action, two-fluid effects introduce dispersion for the ion magnetoacoustic and neutral acoustic waves which are described by Eqs. (1)–(13). As a result of this dispersion, the initial pulse spreads into a train of upward- and downward-propagating waves. The transition region oscillates as a natural consequence of the incoming train of waves which is clearly seen as upward- and downward-propagating waves. These waves move with the local tube speed which is subsonic and subAlfvénic. As gravity introduces characteristic spatial scales in the form of ion and neutral pressure scale heights, the amplitude of the oscillations in V_iy grows with height (top panels). At the top of the photosphere, namely at y  =  0.5 Mm, the signal in V_iy is relatively small and reaches a maximum value there of 0.8 km s⁻¹ for A = 0.5 km s⁻¹ and 9 km s⁻¹ for A = −5 km s⁻¹. Similarly at later times, after t = 10³ s, the signal is very weak as the impulsively generated magnetoacoustic waves have already left the system. Shock waves are generated up to t ≈ 700 s. This means that because of this initial pulse, thermal energy is released until this time because it is mainly generated at the shocks by ion–neutral collisions.

Fig. 3. Time–distance plots for Viy (top) and averaged Viy over time (bottom) for y0 = 0 Mm, A = 0.5 km s−1 (left), and A = −5 km s−1 (right). The velocity is expressed in units of km s−1.

The initial pulse also induces plasma flows. These are evident from the temporarily averaged vertical ion velocity given as

$$\begin{array}{c}\hfill ⟨V_{\mathit{iy}}⟩=\frac{1}{t_2-t_1}{\displaystyle {\int }_{t_1}^{t_2}{V}_{\mathit{iy}}\mathrm{d}t,}\end{array}$$(19)

where t₁ = 0 s and t₂ = 5000 s. The outflows occur in the corona (Fig. 3, bottom). A higher amplitude pulse results in faster outflows originating from lower levels of the atmosphere. In particular, for A = 0.5 km s⁻¹, ⟨V_iy⟩ grows almost linearly with height y, reaching a value of 0.02 km s⁻¹ at y = 20 Mm (Fig. 3, bottom-left). For A = −5 km s⁻¹ the growth of ⟨V_iy⟩ with y is more abrupt than in the case of A = 0.5 km s⁻¹ and ⟨V_iy⟩ attains a value of about 0.7 km s⁻¹ at 20 Mm. This speed is smaller than in the semi-empirical data of Avrett & Loeser (2008), where the velocity is of the order of 40 km s⁻¹ at y = 3 Mm. Downflows are observed in the lower atmospheric regions. Such upflows of ⟨V_iy⟩ = 2 km s⁻¹ and downflows of ⟨V_iy⟩ < 10 km s⁻¹ were reported by Kayshap et al. (2015) and more recently also by Tian et al. (2021).

Figure 4 (left panel) illustrates the ion mass flux, F_𝜚iy = 𝜚_iyV_iy, which in contrast to the ion velocity reaches its highest values below the transition region. Negative values of the flux are observed only in the initial phase of the wave evolution lasting up to t ≈ 200 s. The magnetoacoustic waves, which are reflected from the transition region, are represented by the stripes leaning towards the y-axis. The right panel presents the vertical profile of the averaged (over time) ion mass flux which reaches its lowest value of F_ϱi,y = −0.2 × 10⁻⁶ g cm⁻² s⁻¹ at y = 0.2 Mm. Higher up, it abruptly grows to F_ϱi,y = 10⁻⁶ g cm⁻² s⁻¹ at y = 1.5 Mm and then remains essentially constant with height.

Fig. 4. Time–distance plots for ion mass flux, F𝜚i, y = 𝜚iViy (left), and its profile averaged over time (right) for the initial pulse with width w = 0.3 Mm and amplitude A = 0.5 km s−1, launched from y = y0 = 0 Mm.

Figure 5 (top) illustrates the velocity drift, δV = V_iy − V_ny, which grows with height, reaching its maximum magnitude in the upper chromosphere. However, at the bottom of the numerical region, ions and neutrals are strongly coupled, and therefore they propagate with almost the same velocity. As the quantity δV attains its largest values in the initial phase, we infer that plasma heating can take place up to t = 10³ s after the initial pulse. The values of the velocity drift are small, which indicates that ions and neutrals are strongly coupled. In the right-top panel, at the bottom of the photosphere, y = 0 Mm, and below at y ≈ 0.2 Mm and y ≈ 0.5 Mm we observe some oscillatory signatures. Signatures at y = 0 Mm, and y ≈ 0.2 Mm look like quasi-harmonic decaying oscillations while the third one at y ≈ 0.5 Mm is seen as faster decaying waves propagating upward. As they appear at the bottom of the simulation box they, may result from numerical errors.

Fig. 5. Time–distance plots for the velocity drift, Viy − Vny (top), δTi/T (middle), and averaged δTi/T over time (bottom) for the pulse with width w = 0.3 Mm, launched from y = y0 = 0 Mm and with amplitude: A = 0.5 km s−1 (left) and A = −5 km s−1 (right).

The middle panels of Fig. 5 show the relative perturbed temperature specified as

$$\begin{array}{c}\hfill \frac{\delta T_i}{T}=\frac{T_i-T}{T}.\end{array}$$(20)

This quantity reaches its largest magnitude of 0.55 at y ≈ 2 Mm for A = 0.5 km s⁻¹ and 0.8 at y ≈ 2 Mm for A = −5 km s⁻¹, while the photosphere is hardly heated at all. We note that the signals in the top and middle panels exhibit a level of correlation, which means that the chromospheric plasma is heated by the ion–neutral collisions.

From Eqs. (10) and (11) we infer that the velocity drift, δV, contributes an important part to the heating and this contribution is particularly effective at shocks, which is clearly seen in Fig. 5 (top). This was also noticed by Zhang et al. (2021) and Snow & Hillier (2021) from similar simulations in a different setup. This heating can be quantified by the quantity δT_i/T which is displayed on Fig. 5 (middle).

We define the time-averaged relative perturbed temperature as

$$\begin{array}{c}\hfill \langle \frac{\delta T_i}{T}\rangle =\frac{1}{t_2-t_1}{\displaystyle {\int }_{t_1}^{t_2}\frac{\delta T_i}{T}dt,}\end{array}$$(21)

where t₁ = 0 s and t₂ = 5000 s. This quantity is illustrated in the bottom panel of Fig. 5. We note that the chromosphere is more strongly heated for higher values of the pulse amplitude |A|. The lowest averaged temperature is seen in the top of the photosphere and above this region the atmosphere is heated to a greater temperature. For A = −5 km s⁻¹, all values of ⟨δT_i/T⟩ are positive and the minimum of ⟨δT_i/T⟩ ≈ 0 at y ≈ 0.4 Mm takes place at the equilibrium temperature minimum (Avrett & Loeser 2008), where the magnetoacoustic waves experience reflections. ⟨δT_i/T⟩ reaches a value of about 2.3 at y = 2 Mm, at the top of the chromosphere. A characteristic flattening is present within the lower levels of the middle chromosphere, from y ≈ 0.8 Mm to y ≈ 1 Mm. Above this plateau, ⟨δT_i/T⟩ suddenly increases with y (bottom-right).

We note that the temperature is reduced permanently in the photosphere with essentially unseen waves in δT_i/T (Fig. 5, middle). In the chromosphere, we observe oblique stripes which are leaned towards the vertical axis and represent magnetoacoustic waves. In the case of A = 0.5 km s⁻¹, max (δT_i/T) ≈ 0.6 at y = 1.8 Mm, which is located in the chromosphere (middle-left). As for larger absolute values of amplitude, for example for A = −5 km s⁻¹ (middle-right), the generated shocks are more pronounced than for A = 0.5 km s⁻¹, and result in more intensive heating with max (δT_i/T) ≈ 0.8 in the middle of the chromosphere at y = 1.1 Mm (Fig. 5, middle).

3.2.2. Pulse width effect

In order to study how the width of the pulse, w, affects the plasma heating, we performed some parametric studies. From Fig. 6 (top) we infer that increasing w from 0.1 Mm to 0.25 Mm results in a significant heating of the chromosphere. However, for w = 0.1 Mm, only the top layer of the chromosphere is heated (top-left), while a wider pulse width, w = 0.3 Mm, results in heating of the whole chromosphere (top-right). Similarly, the cooling during the first 200 s is more pronounced for the wider pulse, and results from the plasma being rarefied.

Fig. 6. Time–distance plots for δTi/T for the pulse with amplitude A = 2 km s−1 and width w = 0.1 Mm (top-left) and w = 0.25 Mm launched from y = y0 = 0 Mm (top-right), and for w = 0.3 Mm launched from y0 = 0.25 Mm (bottom-left) and y0 = 0.5 Mm (bottom-right).

3.2.3. Pulse launching height effect

Figure 6 (bottom) shows time–distance plots for δT_i/T for different pulse-launching heights, that is, for different values of the parameter y₀. A higher value of y₀ results in a smaller value of max (δT_i/T) at a given height y. This is the consequence of a pressure scale-height over which, according to linear theory, the amplitude of the signal grows e−times (e.g., Nakariakov & Verwichte 2005). For the pulse launched from a higher level y₀, with the same amplitude, the signal passes fewer scale heights, reaching a given level y. Consequently, δV and δT_i/T are smaller (at a given height y) for a larger value of y₀.

Figure 7 illustrates the relative perturbed temperature averaged over time and height, H, defined as

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

Fig. 7. Relative perturbed temperature averaged over time and height, H, vs. |A|, for y0 = 0 Mm and w = 0.3 Mm (top); w for A = 2 km s−1 and y0 = 0 Mm (middle); and y0 for A = 2 km s−1 and w = 0.3 Mm (bottom).

Here, y₁ = 1.9 Mm.

The higher value of the pulse amplitude |A| results in more significant heating of the photosphere and chromosphere (top). Similarly, a higher magnitude of the pulse width w leads to a higher value of H (middle). Finally, in creasing the pulse launching height y₀ corresponds to less heating below the transition region with a fall-off of H(y₀) resulting from the pressure scale height (bottom). For a larger value of y₀, the distance to the transition region covers fewer pressure scale heights and the signal has less chance to increase its amplitude and consequently to generate stronger shocks with less thermal energy release. The value of H = 1 would correspond to an average relative temperature increase in the photosphere and chromosphere by 100%. However, we do not reach such a value in our results because the amplitude of the pulse is insufficient to induce such strong heating.

3.2.4. Outflow velocities and Fourier power spectrum

Figure 8 shows the maximum ion outflow velocity, max(V_iy), for different pulse parameters. Similarly to in Fig. 7, a pulse with larger absolute amplitude (top) and width (middle) results in a higher value of the maximum outflow velocity max(V_iy). As we increase the launching altitude of the initial pulse, y₀, the velocity naturally decreases (bottom).

Fig. 8. Maximum value of Viy vs. |A|, for y0 = 0 Mm and w = 0.3 Mm (top); w for A = 2 km s−1 and y0 = 0 Mm (middle); and y0 for A = 2 km s−1 and w = 0.3 Mm (bottom).

Figure 9 presents Fourier power spectrum for wave period P versus height. The initial pulse has a Gaussian spectrum of wave number k which results in a spectrum of period P. The steepening of the magnetoacoustic waves results from the growing wave amplitude with height. Hence, waves with shorter wavelengths and wave periods are present for higher y-values in their Fourier spectra. For the pulse launched from y = y₀ = 0 Mm (top), the main wave period of the downward-propagating waves becomes approximately the same for all values of y < 0 Mm and is equal to about 250 s. Higher up, however, the period P decays with increasing y, and attains values close to 200 s. As a result of the cut-off only short-period waves can propagate upwards while long-period waves become evanescent. Hence, the relative contribution of long P waves weakens with increasing y. We note that some of our data fits the observational findings of Wiśniewska et al. (2016), represented by diamonds over-plotted on the power spectra, and Kayshap et al. (2018), denoted by dots. The agreement of the theory with the observational data indicates that the results can be used to determine the background structure of the solar atmosphere and confirms that wave generation by the solar granulation in the partially ionized plasma dominates the behavior of the waves. For a pulse at the bottom of the photosphere (y = y₀ = 0 Mm, Fig. 9 (top)), a jump in the dominant wave period is observed close to the height y = 1.5 Mm, which reaches a magnitude of 300 s. The signal launched from the bottom of the photosphere with the main period P = 250 s thus reaches the corona with P = 300 s. The wave periods in the photosphere are lower than P_ac in this region, because max P_ac = 240 s at y = 0.5 Mm, which means that magnetoacoustic waves are evanescent. However, in the corona above y = 2.3 Mm the acoustic cut-off period reaches values larger than P = 300 s. However, this is not the case for a pulse launched at a somewhat greater height, in the middle of the photosphere, namely for y = y₀ = 0.25 Mm (see Fig. 9 (bottom)). Moreover, in this case the layers below the photosphere (y < 0) oscillate with a dominant wave period of 225 s. In the photosphere, this wave period is also dominant and this time it increases slightly with height in the photosphere. Moreover, in this case there is a second dominant period, namely 250 s. In fact, this wave period is also the dominant one in the upper chromosphere and low corona in this case.

Fig. 9. Fourier period P for Viy vs. height, in the case of the initial pulse with amplitude A = 0.5 km s−1, width w = 0.3 Mm, and launched from y = y0 = 0 Mm (top) and y = y0 = 0.25 Mm (bottom). 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. Summary and conclusions

In this paper we consider two-fluid ion magnetoacoustic and neutral acoustic waves that are initially (at >t = 0 s) excited by a pulse in a vertical component of the ion and neutral velocities. The full set of two-fluid equations is solved by the JOANNA code (Wójcik et al. 2020). The results of the numerical simulations presented in this paper can be summarized as follows.

The triggered magnetoacoustic waves transform into shock waves and through ion–neutral collisions they convert their energies into thermal energy contributing to heating of the chromospheric plasma. Increasing the amplitude of the pulse leads to a significant heating of the chromosphere, which for the highest absolute value of the amplitude considered, A = −5 km s⁻¹, is clearly seen to occur in the upper part of the chromosphere, above the level y = 1 Mm. Moreover, the heating of the chromosphere is correlated with a velocity drift which increases with height up to y = 2 Mm. This means that the chromosphere is heated by the ion–neutral collisions. The parametric studies we performed for various widths of the pulse show that a wider pulse corresponds to more plasma heating spread over the whole chromosphere. On the other hand, the investigation of quasi-periodic MHD waves by Nakariakov et al. (2005) showed that a wider pulse is less effective in generating quasi-periodic wave trains. The problem of impulsively generated wave trains was also studied by Yu et al. (2017) who concluded that the diversity of group speed characteristics has an impact on the temporal evolution of impulsively generated wave trains. A pulse launched from a higher level (with the same amplitude), on the other hand, results in less heating as the amplitude of upwardly propagating signal cannot grow significantly with height before reaching the transition region. Both the plasma temperature and maximum ion velocity increase with increasing pulse amplitude and width, while the opposite trend is seen for higher values of y₀, that is, for launching the pulse higher up. The pulses considered in this paper result in shocks that lead to a release of thermal energy which in turn helps to sustain the chromosphere. We arrive at a further conclusion, namely that even a single pulse can drive vertical plasma outflows which higher up can result in (contributions to) the origin of the solar wind.

Acknowledgments

The authors express their gratitude to the reviewer for his/her stimulating comments and suggestions. The JOANNA code was written by Dr. Dariusz Wójcik. This work was done within the framework of the project from the National Science Centre (NCN) grant nos. 2017/25/B/ST9/00506 and 2020/37/B/ST9/00184. Numerical simulations were performed on the LUNAR cluster at the Institute of Mathematics of University of M. Curie-Skłodowska, Lublin, Poland. SP is supported by the project EUHFORIA 2.0 which has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 870405, and by the projects C14/19/089 (C1 project Internal Funds KU Leuven), G.0D07.19N (FWO-Vlaanderen), SIDC Data Exploitation (ESA Prodex-12), and the Belspo projects BR/165/A2/CCSOM and B2/191/P1/SWiM.

References

All Figures

	Fig. 1. Variation with height of the equilibrium temperature (top), mass densities (middle) of neutrals (dashed line) and ions (dotted line), and acoustic cut-off period (bottom).
In the text

	Fig. 2. Time–distance plots for Viy/Cs (top) and δTi/T (bottom) in the case of no pulse (A = 0).
In the text

	Fig. 3. Time–distance plots for Viy (top) and averaged Viy over time (bottom) for y0 = 0 Mm, A = 0.5 km s−1 (left), and A = −5 km s−1 (right). The velocity is expressed in units of km s−1.
In the text

	Fig. 4. Time–distance plots for ion mass flux, F𝜚i, y = 𝜚iViy (left), and its profile averaged over time (right) for the initial pulse with width w = 0.3 Mm and amplitude A = 0.5 km s−1, launched from y = y0 = 0 Mm.
In the text

	Fig. 5. Time–distance plots for the velocity drift, Viy − Vny (top), δTi/T (middle), and averaged δTi/T over time (bottom) for the pulse with width w = 0.3 Mm, launched from y = y0 = 0 Mm and with amplitude: A = 0.5 km s−1 (left) and A = −5 km s−1 (right).
In the text

	Fig. 6. Time–distance plots for δTi/T for the pulse with amplitude A = 2 km s−1 and width w = 0.1 Mm (top-left) and w = 0.25 Mm launched from y = y0 = 0 Mm (top-right), and for w = 0.3 Mm launched from y0 = 0.25 Mm (bottom-left) and y0 = 0.5 Mm (bottom-right).
In the text

	Fig. 7. Relative perturbed temperature averaged over time and height, H, vs. |A|, for y0 = 0 Mm and w = 0.3 Mm (top); w for A = 2 km s−1 and y0 = 0 Mm (middle); and y0 for A = 2 km s−1 and w = 0.3 Mm (bottom).
In the text

	Fig. 8. Maximum value of Viy vs. |A|, for y0 = 0 Mm and w = 0.3 Mm (top); w for A = 2 km s−1 and y0 = 0 Mm (middle); and y0 for A = 2 km s−1 and w = 0.3 Mm (bottom).
In the text

	Fig. 9. Fourier period P for Viy vs. height, in the case of the initial pulse with amplitude A = 0.5 km s−1, width w = 0.3 Mm, and launched from y = y0 = 0 Mm (top) and y = y0 = 0.25 Mm (bottom). 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