EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 691 (November 2024) A&A, 691 (2024) A254 Full HTML
Open Access
Issue
A&A
Volume 691, November 2024
Article Number A254
Number of page(s) 7
Section The Sun and the Heliosphere
DOI https://doi.org/10.1051/0004-6361/202449941
Published online 18 November 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

The Sun’s atmosphere is classified into three layers that differ in physical quantities, such as temperature and mass density. The lowest layer of the atmosphere is called the photosphere, where the temperature ranges from about 5600 K at its bottom to about 4300 K at its top, located about 600 km above the solar surface. Above the photosphere is the chromosphere, and in this layer, the temperature rises to about 6 − 7 ⋅ 103 K. In the outer layer, called the solar corona, the temperature reaches 1 − 3 MK on average. As the abundance of ions is strongly correlated with the plasma temperature, the photosphere is only very weakly ionised, that is, the ionisation degree is of the order of 10−4. The chromosphere, however, is partially ionised with typical ionisation levels of 5 ⋅ 10−2, and the corona is essentially fully ionised (Biermann 1947; Miyamoto 1949). The solar corona and the chromosphere are separated by a thin layer (100 − 200 km) called the transition region, in which temperature suddenly rises from chromospheric to coronal values. This abundance of neutrals in the lower atmospheric layers means that the presence of both charged species (ions and electrons) and neutral atoms must be considered to describe the dynamics of the atmosphere realistically.

Being a dynamic environment, the photosphere is the source of diverse waves. The effect of these waves on the chromospheric temperature and the generation of plasma flows has been investigated in recent years (e.g. Nakariakov et al. 2019; Niedziela et al. 2022; Pelekhata et al. 2022). More realistic models of the atmosphere, including complex magnetic fields and non-isothermal conditions, have been considered in a series of studies (e.g. Defouw 1976; Stark & Musielak 1993; Musielak et al. 2006; Felipe et al. 2018; Murawski et al. 2022; Kraśkiewicz et al. 2022).

Naturally, ion–neutral collisions result from the chromosphere’s physical conditions. A two-fluid model could model them, but the effect of these collisions can be partly taken into account in a magnetohydrodynamic (MHD) approach, such as ambipolar diffusion in the generalised Ohm’s law (e.g. MacBride et al. 2022; González-Morales et al. 2020). A full two-fluid model is superior to an MHD model, and it has been shown that these ion–neutral collisions may play a key role in wave damping – a process that is most effective at steep wave profiles (e.g. Kuźma et al. 2021; Murawski et al. 2020).

Driven slow magnetoacoustic waves (SMAWs) were also considered in the framework of MHD by Kraśkiewicz et al. (2019), who performed numerical simulations of the behaviour of excited SMAWs in MHD systems with horizontal and vertical magnetic fields. Further numerical studies of SMAWs generated by the solar granulation show that the main period of highly turbulent photospheric plasma is close to 300 s (Kuźma et al. 2021).

The idea of studying chromosphere heating by acoustic waves was initially suggested by Biermann (1946) and Schwarzschild (1948). Much later, it was proposed that the incompressible counterparts of acoustic waves, mainly Alfvén waves, can participate in the generation of the solar wind (e.g. Ofman 2010; Banerjee et al. 2021).

The main aim of the present paper is to study the propagation of driven two-fluid SMAWs in the partially ionised solar atmosphere. The solar atmosphere models developed thus far (e.g. Botha et al. 2011; Kraśkiewicz et al. 2023; Niedziela et al. 2021; Pelekhata et al. 2021) have been based on a pure hydrostatic equilibrium (HE) approximation, whereas the model presented here supplements this for the first time with the Saha equation (SE) (Saha 1920) in order to consider the ionisation effects. This equation is important as it improves the HE model by implementing a realistic vertical profile of the ionisation degree, that is, the ratio of the number of ions to the number of neutrals in the solar atmosphere.

This paper is organised as follows. The numerical model and the setup of the numerical simulations are described in Sect. 2. The numerical results of the driven waves are presented and discussed in Sect. 3. Our findings and conclusions are summarised in Sect. 4.

2. Numerical model

To describe the lower layers of the solar atmosphere, we considered a hydrogen plasma, which consists of two components: an ionised (protons + electrons) fluid and a neutral (hydrogen atoms) fluid.

2.1. Two-fluid equations

The appropriate two-fluid equations are given by Oliver et al. (2016), Maneva et al. (2017), Popescu Braileanu et al. (2019), Murawski et al. (2022) and rewritten here as

∂ ϱ i ∂ t + ∇ · ( ϱ i V i ) = m i ( Γ i ion + Γ i rec ) , $$ \begin{aligned}&\frac{\partial \varrho _{\rm i}}{\partial t}+\nabla \cdot (\varrho _{\rm i} \boldsymbol{V}_{\rm i}) = m_{\rm i} (\Gamma ^\mathrm{ion}_{\rm i}+\Gamma ^\mathrm{rec}_{\rm i}) , \end{aligned} $$(1)

∂ ϱ n ∂ t + ∇ · ( ϱ n V n ) = m n ( Γ n ion + Γ n rec ) , $$ \begin{aligned}&\frac{\partial \varrho _{\rm n}}{\partial t}+ \nabla \cdot (\varrho _{\rm n} \boldsymbol{V}_{\rm n}) = m_{\rm n} (\Gamma ^\mathrm{ion}_{\rm n}+\Gamma ^\mathrm{rec}_{\rm n}) , \end{aligned} $$(2)

∂ ( ϱ i V i ) ∂ t + ∇ · ( ϱ i V i V i + ( p i + p e ) I ) = ϱ i g + 1 μ ( ∇ × B ) × B + S m , $$ \begin{aligned}&\frac{\partial (\varrho _{\rm i} \boldsymbol{V}_{\rm i})}{\partial t} + \nabla \cdot (\varrho _{\rm i} \boldsymbol{V}_{\rm i} \boldsymbol{V}_{\rm i} + ( p_{\rm i} +p_{\rm e}) \boldsymbol{I}) = \varrho _{\rm i} \boldsymbol{g} + \frac{1}{\mu }(\nabla \times \boldsymbol{B}) \times \boldsymbol{B} + \boldsymbol{S_{\rm m}} , \end{aligned} $$(3)

∂ ( ϱ n V n ) ∂ t + ∇ · ( ϱ n V n V n + p n I ) = ϱ n g − S m , $$ \begin{aligned}&\frac{\partial (\varrho _{\rm n} \boldsymbol{V}_{\rm n})}{\partial t} + \nabla \cdot (\varrho _{\rm n} \boldsymbol{V}_{\rm n} \boldsymbol{V}_{\rm n} + p_{\rm n} \mathbf I ) = \varrho _{\rm n} \boldsymbol{g} - \boldsymbol{S_{\rm m}} , \end{aligned} $$(4)

∂ E i ∂ t + ∇ · [ ( E i + p i + p e + B 2 2 μ ) V i − B μ ( V i · B ) ] = ( ϱ i g + S m ) · V i + Q i , $$ \begin{aligned}&\frac{\partial E_{\rm i}}{\partial t}+\nabla \cdot \left[\left(E_{\rm i}+p_{\rm i} + p_{\rm e} + \frac{{\boldsymbol{B}}^2}{2\mu } \right)\boldsymbol{V}_{\rm i}-\frac{{\boldsymbol{B}}}{\mu }(\boldsymbol{V}_{\rm i}\cdot \boldsymbol{B}) \right] = (\varrho _{\rm i} \boldsymbol{g}+\boldsymbol{S_{\rm m}}) \cdot \boldsymbol{V}_{\rm i} + Q_{\rm i}, \end{aligned} $$(5)

∂ E n ∂ t + ∇ · [ ( E n + p n ) V n ] = ( ϱ n g − S m ) · V n + Q n , $$ \begin{aligned}&\frac{\partial E_{\rm n}}{\partial t}+\nabla \cdot [(E_{\rm n}+p_{\rm n})\boldsymbol{V}_{\rm n}] = (\varrho _{\rm n} \boldsymbol{g} - \boldsymbol{S_{\rm m}}) \cdot \boldsymbol{V}_{\rm n} + Q_{\rm n} , \end{aligned} $$(6)

∂ B ∂ t = ∇ × ( V i × B ) , ∇ · B = 0 , $$ \begin{aligned}&\frac{\partial \boldsymbol{B}}{\partial t} = \nabla \times (\boldsymbol{V_{\rm i} \times }\boldsymbol{B}) , \quad \nabla \cdot \boldsymbol{B} = 0 , \end{aligned} $$(7)

E i = ϱ i V i 2 2 + p i + p e γ − 1 + B 2 2 μ , and E n = ϱ n V n 2 2 + p n γ − 1 . $$ \begin{aligned}&E_{\rm i} = \frac{\varrho _{\rm i}\boldsymbol{V}_{\rm i}^2}{2} + \frac{p_{\rm i}+p_{\rm e}}{\gamma -1 } +\frac{\boldsymbol{B}^2}{2\mu } , \quad \mathrm{and} \quad E_{\rm n} = \frac{\varrho _{\rm n}\boldsymbol{V}_{\rm n}^2}{2} + \frac{p_{\rm n}}{\gamma -1 } . \end{aligned} $$(8)

In the equations above, the subscripts n, i, and e correspond to neutrals, ions, and electrons, respectively, while 𝜚i, n denotes the mass densities of ions and neutrals, Vi, n are the velocities, and pi, e, n describe the ion, electron, and neutral gas pressures. It is assumed that pe = pi. The symbol I stands for the identity matrix, the terms Ei, n correspond to the total energy densities, g = [0, −g, 0] is the gravity vector with magnitude g = 274.78 m s−2, γ = 5/3 represents the ratio of specific heats, B denotes the magnetic field, and μ is the magnetic permeability. The symbols Sm and Qi, n indicate the collisional momentum and the energy exchange terms. These are given by the following expressions (Oliver et al. 2016):

S m = α in ( V n − V i ) + Γ i ion m i V n − Γ n rec m i V i , $$ \begin{aligned} \boldsymbol{S}_{\rm m}&= \alpha _{\rm in} (\boldsymbol{V_{\rm n}} - \boldsymbol{V}_{\rm i}) + \Gamma ^\mathrm{ion}_{\rm i} m_{\rm i} \boldsymbol{V}_{\rm n} - \Gamma ^\mathrm{rec}_{\rm n} m_{\rm i} \boldsymbol{V}_{\rm i} , \end{aligned} $$(9)

Q i = 1 2 α in ( V i − V n ) 2 + 1 γ − 1 k B α in m n + m i ( T n − T i ) + 1 2 m i ( Γ i ion V n 2 − Γ i rec V i 2 ) , and $$ \begin{aligned} Q_{\rm i}&= \frac{1}{2}\alpha _{\rm in} (\boldsymbol{V}_{\rm i}-\boldsymbol{V}_{\rm n})^2 + \frac{1}{\gamma -1} \frac{k_{\rm B}\alpha _{\rm in}}{m_{\rm n}+m_{\rm i}} (T_{\rm n}-T_{\rm i}) \nonumber \\&+\frac{1}{2} m_{\rm i} (\Gamma ^\mathrm{ion}_{\rm i} V^2_{\rm n} - \Gamma ^\mathrm{rec}_{\rm i} V^2_{\rm i}) ,\,\, \mathrm{and} \end{aligned} $$(10)

Q n = 1 2 α in ( V i − V n ) 2 + 1 γ − 1 k B α in m n + m i ( T i − T n ) + 1 2 m i ( Γ n rec V i 2 − Γ n ion V n 2 ) . $$ \begin{aligned} Q_{\rm n}&= \frac{1}{2}\alpha _{\rm in} (\boldsymbol{V_{\rm i}}-\boldsymbol{V_{\rm n}})^2 + \frac{1}{\gamma -1}\frac{k_{\rm B}\alpha _{\rm in}}{m_{\rm n}+m_{\rm i}} (T_{\rm i}-T_{\rm n}) \nonumber \\&+\frac{1}{2} m_{\rm i} (\Gamma ^\mathrm{rec}_{\rm n} V^2_{\rm i} - \Gamma ^\mathrm{ion}_{\rm n} V^2_{\rm n}) . \end{aligned} $$(11)

The ionisation and recombination coefficients, Γion and Γrec, are specified by Maneva et al. (2017). The effect of the interaction between these two species depends on the ion–neutral friction coefficient αin, which is defined as (Braginskii 1965)

α in = 4 3 σ in ϱ i ϱ n m i + m n 8 k B π ( T i m i + T n m n ) . $$ \begin{aligned} \alpha _{\rm in} = \frac{4}{3} \frac{\sigma _{\rm in}\varrho _{\rm i}\varrho _{\rm n}}{m_{\rm i}+m_{\rm n}} \sqrt{ \frac{8k_{\rm B}}{\pi } \left( \frac{T_{\rm i}}{m_{\rm i}}+\frac{T_{\rm n}}{m_{\rm n}} \right)} . \end{aligned} $$(12)

Here, kB corresponds to the Boltzmann constant, σin = 1.4 ⋅ 10−15 cm2 represents the quantum collision cross-section (Vranjes & Krstic 2013), mi, n are the masses of the ions (protons) and the neutrals (hydrogen atoms), respectively, and Ti, n are their temperatures given by the ideal gas laws,

p i , n = k B m i , n ϱ i , n T i , n . $$ \begin{aligned} p_{\rm i,n}= \frac{k_{\rm B}}{m_{\rm i,n}}\varrho _{\rm i,n}T_{\rm i,n}. \end{aligned} $$(13)

For the sake of simplicity, all other non-adiabatic and non-ideal terms in the two-fluid equations were neglected, and the discussion is limited to a 2D situation with z being an invariant coordinate. However, compressive viscosity, which was not included in the model, would additionally damp the magnetoacoustic waves (Nakariakov et al. 2017; Duckenfield et al. 2021; Ofman & Wang 2022).

2.2. Magnetohydrostatic equilibria

Two realisations of the solar atmosphere equilibrium were considered, namely:

  • 1.

    a HE, which was determined by a vertical temperature profile T0(y), taken here from the model of Avrett & Loeser (2008), with reference values of ion and neutral gas pressures (for details, see e.g. Niedziela et al. 2022);

  • 2.

    a HE supplemented by the SE (HE+SE). Here, the hydrostatic equations were combined with the SE, which is specified for the hydrogen plasma as (Suzuki et al. 2022)

    n i + 1 n i = 1 n H λ e 3 exp ( − I H k B T ) . $$ \begin{aligned} \frac{n_{\rm i+1}}{n_{\rm i}} = \frac{1}{n_{\rm H} \lambda _{\rm e}^3} \exp \Big (-\frac{I_{\rm H}}{k_{\rm B} T}\Big ) . \end{aligned} $$(14)

    Here, ni + 1 and ni are the density number of atoms in the ionisation state i + 1 and i, respectively, nH denotes the density of hydrogen atoms, which in the hydrogen plasma model is equal to nn, IH = 13.6 eV corresponds to the hydrogen ionisation potential, and λe is the thermal de Broglie wavelength of an electron, given as

    λ e = h 2 2 π m e k B T . $$ \begin{aligned} \lambda _{\rm e} = \sqrt{\frac{h^2}{2 \pi m_{\rm e} k_{\rm B} T}} . \end{aligned} $$(15)

    Here, h represents the Planck constant, and me is the electron mass. Note that the SE is valid for a plasma in thermodynamic equilibrium, which is not valid for the solar atmosphere.

The vertical profile of the equilibrium temperature, T0(y), is displayed up to y = 5 Mm in Fig. 1a. In the photosphere, T0 attains a value of about 5600 K. At y = 0.6 Mm, which is close to the bottom of the chromosphere, it decreases to a minimum of about 4300 K (Fig. 1a). In contrast, T0(y) grows slightly in the chromosphere, but the most dramatic variation is observed in the transition region, across which the temperature rises from about 7 ⋅ 103 K to 2 ⋅ 105 K at y = 2.5 Mm.

The ionisation degree, I = 𝜚i/(𝜚i + 𝜚n), determined with the use of the SE reached the lowest values at the temperature minimum at the top of the photosphere, with a minimum value of about 2 ⋅ 10−5 (Fig. 1, panel b, solid line). Moreover, it was significantly different for the HE case, for which I attained the lowest value of about 10−2 at y = 0 Mm (dashed line). As T0(y) increases in the chromosphere, I also tends to grow there. Due to the essentially complete ionisation of plasma in the solar corona resulting from its high temperature, the magnitudes of I are highest there. The maximum value of I ≈ 1 takes place at y = 2.3 Mm and corresponds to the floor value for 𝜚n, which was set as 10−19 g ⋅ cm−3 to avoid numerical issues with negative mass density in the corona. Note that 𝜚i is lower than 𝜚n in the photosphere and the chromosphere (Fig. 1c). This results from the relatively low chromosphere temperature compared to the temperature of the solar corona, which directly affects the ionisation rate and the number of ions in the system. At y = 0.5 Mm and at y ≈ 1 Mm, respectively, a local minimum and maximum of 𝜚i take place in the HE+SE case. However, in the transition region, where a sudden increase in temperature occurs, a sharp decrease in 𝜚n is observed. In the corona, the mass density reaches its minimum, which is higher for ions, with 𝜚i ≈ 10−14 g cm−3, than for neutrals, with 𝜚n ≈ 10−19 g cm−3 for HE and 𝜚n ≈ 10−16 g cm−3 for HE+SE (Figs. 1c,d). We clearly see that for HE, 𝜚i already dominates over 𝜚n in the middle chromosphere, that is, at a lower level than for HE+SE. The hydrostatic system is embedded in a magnetic field,

thumbnail Fig. 1.

Variation with height of the equilibrium temperature, T0 (a), ionisation degree, I (b) for HE (dashed line) and HE+SE (solid line), and mass density (c, d) for ions (solid line) and neutrals (dashed line) for HE+SE (c) and HE (d).

B = B 0 [ cos ( x + L Λ B ) , − sin ( x + L Λ B ) , 0 ] exp ( − y Λ B ) + [ 0 , B v , 0 ] , $$ \begin{aligned} \boldsymbol{B}&= B_{\rm 0}\left[ \cos \left(\frac{x+L}{\Lambda _{\rm B}} \right), -\sin \left( \frac{x+L}{\Lambda _{\rm B}} \right),0\right] \exp \left( -\frac{{ y}}{\Lambda _{\rm B}} \right) \nonumber \\ &+ [0,B_{\rm v},0],\end{aligned} $$(16)

Λ B = 2 L / π , $$ \begin{aligned} \Lambda _{B}&= 2 L/\pi , \end{aligned} $$(17)

which mimics fanning out magnetic field lines (Fig. 2). The magnitude Bv of the vertical field is chosen as −5 Gs. The magnetic field at y = 0, B0 = 500 Gs, and the characteristic width L = 0.64 Mm. This magnetic field is current-free (∇×B/μ = 0) and thus also force-free ((∇ × B) × B/μ = 0), and therefore it does not affect the hydrostatic ion gas pressure and ion mass density profiles.

thumbnail Fig. 2.

Spatial profile of Vi(x, y) (expressed in km s−1) at t = 110 s overlaid by its vectors (white arrows) and magnetic field lines for Pd = 5 s.

2.3. Perturbations by a monochromatic driver

The magnetohydrostatic equilibrium is perturbed by the monochromatic driver in vertical components of ion and neutral velocities, given as

V i y ( y = y d , t ) = V i y ( y = y d , t ) = A exp ( − x 2 w 2 ) sin ( 2 π t P d ) . $$ \begin{aligned} V_{\mathrm{i}{ y}}({ y}={ y}_{\rm d},t)=V_{\mathrm{i}{ y}}({ y}={ y}_{\rm d},t)= A \exp {\left(-\frac{x^2}{{ w}^2} \right)} \sin \left(\frac{2 \pi t}{P_{\rm d}}\right). \end{aligned} $$(18)

Here, A is the amplitude of the driver, which is fixed to 5 km s−1, w = 50 km its width, Pd is its period, and yd denotes the operation height which is set to be equal to yd = 1 Mm, corresponding to the middle chromosphere. These amplitudes were observed in the chromosphere in three-minute oscillations (Felipe et al. 2010; Krishna Prasad et al. 2015; Khomenko & Collados 2015). Several values of Pd were considered, and appropriate results for the HE were compared with those for the HE+SE model. The following values of Pd = 2.5 s, 5 s, and 10 s were considered. As these values are closer to the ion–neutral collision times, which are within the range of about 0.1 − 1 s in the chromosphere and the low corona (Khomenko et al. 2014), the two-fluid effects were anticipated to be greater for these periods than for the more powerful three- and five-minute oscillations.

2.4. Numerical box and boundary conditions

The numerical simulations were performed using the JOANNA code (Wójcik et al. 2020), which solves the two-fluid Equations (1)–(13). The Courant-Friedrichs-Levy (CFL) number was set equal to 0.9, and the third-order strong stability preserving Runge-Kutta method (Durran 2010) was adopted and supplemented by the Harten-Lax-van Leer Discontinuities (HLLD) approximate Riemann solver (Miyoshi & Kusano 2005). The 2D simulation domain was specified as (−0.64 ≤ x ≤ 0.64) Mm × (1 ≤ y ≤ 3.56) Mm and covered by 1024 × 2048 cells, leading to a numerical grid cell size in this area Δy = 0.125 km. Higher up, within the zone 3.56 Mm ≤ y ≤ 15 Mm, the grid was stretched along y up to y = 15 Mm and covered by 64 cells. This stretched grid damped the incoming signal at the top boundary (Kuźma & Murawski 2018). The selection of the numerical grid was preceded by a grid convergence study. The assumption that at least 16 grid points should cover the wavelength proved adequate, leading to numerical diffusion, which did not significantly affect the results for the chosen periods. All plasma quantities were fixed to their equilibrium values at the top and bottom boundaries of the numerical box. The only exception was the bottom boundary overlaid by the monochromatic driver, described by Eq. (18). At the side boundaries, x = ±0.64 Mm, open boundary conditions were implemented by copying all two-fluid quantities into the boundary cells from the nearest physical cells.

3. Numerical results

In this part of the paper, we investigate various driver periods and how they influence the chromosphere heating and the generation of plasma outflows.

Figure 2 illustrates the Vi(x, y) profile at t = 110 s for Pd = 5 s. We note that the leading signal has already passed the transition region and reached y ≈ 2.3 Mm. In the middle and top chromosphere, the periodic pattern can be clearly seen. The signal in Vi is essentially associated with SMAWs, which propagate mainly along magnetic field lines in a strongly magnetised medium.

Figure 3 displays time–distance plots for Viy (a, c, e) and

δ T i / T 0 = T i − T 0 T 0 $$ \begin{aligned} \delta T_{\rm i}/T_{\rm 0}= \frac{T_{\rm i}-T_{\rm 0}}{T_{\rm 0}} \end{aligned} $$(19)

(b, d, f) for Pd = 10 s (a, b), Pd = 5 s (c, d), and Pd = 2.5 s (e, f), collected at x = 0. As a result of the ion–neutral collisions, a part of the SMAWs kinetic energy is thermalised, and this effect is particularly important at steep wave profiles, which result from the exponentially growing wave amplitudes over a height equal to the pressure-scale height. A competitive effect of this is wave energy spreading along magnetic field lines that fan out at height. Thus, even a small amplitude driver may result in SMAWs with quickly growing amplitudes that steepen into shock waves in the chromosphere and thermalise their energy there. The case of Pd = 2.5 s corresponds to the lowest values of max(Viy) = 3.2 km s−1. However, this case shows the greatest heating in the lower part of the numerical domain, below y = 1.2 Mm. The driver with a period of Pd = 5 s excites higher Viy values but lower δTi/T0 compared to the Pd = 2.5 s case. The largest velocities are observed for Pd = 10 s driver. Nevertheless, the lowest heating corresponds to the largest period studied in this paper. Short period waves (Pd = 2.5 s) experience stronger non-linear damping (Kraśkiewicz et al. 2023), which may be an explanation for the lower max(Viy). The incoming signal is partially reflected from the transition region. This reflection results from the physical nature of this layer, where mass density, gas pressure, and temperature experience a sudden fall-off with height.

thumbnail Fig. 3.

Time-distance plots for Viy(x = 0, y, t) (a, c, e) and δTi/T0 (x = 0, y, t) (b, d, f) for Pd = 10 s (a, b), Pd = 5 s (c, d), and Pd = 2.5 s (e, f), in the case of HE+SE.

Figure 4 shows the kinetic energy flux FEics = 0.5𝜚V2CS (a, c) and the frictional heating term Q (b, d), defined as the first term on the right-hand side of Eq. (10). The SMAWs can carry kinetic energy with the flux of about 6 ⋅ 103 erg cm−2 s−1 in the middle chromosphere, where the driver operates. This value falls off with height reaching about 102 erg cm−2 s−1 in the top chromosphere. Compared to the chromospheric radiative energy losses for the quiet Sun (Withbroe & Noyes 1977), the kinetic energy flux is too small to fulfil these losses by about two orders of magnitude. As Q attains the largest values of Q = 2 ⋅ 10−5 erg cm−3 s−1 above the driver, at y ≈ 1.05 Mm (b), we infer that the thermal energy release takes place at this altitude in the considered model atmosphere with active ionisation/recombination effects.

thumbnail Fig. 4.

Time–distance plots (a, b) and averaged-over-time vertical profiles (c, d) for the kinetic energy flux FEics (a, c) and frictional heating term Q (b, d) for Pd = 5  s, in the case of HE+SE.

Figure 5 presents the temporal averaged δTi/T0 (a, c, e) and Viy (b, d, f), defined as

〈 f 〉 t = 1 t 2 − t 1 ∫ t 1 t 2 f d t , $$ \begin{aligned} \Big \langle f \Big \rangle _t = \frac{1}{t_{\rm 2} - t_{\rm 1}} \int ^{t_{\rm 2}}_{t_{\rm 1}} f \, \mathrm{d}t, \end{aligned} $$(20)

thumbnail Fig. 5.

Temporarily averaged ⟨δTi/T0⟩t (a, c, e) and ⟨Viy⟩t (b, d, f) for Pd = 10 s (a, b), Pd = 5 s (c, d), and Pd = 2.5 s (e, f), in the case of HE+SE.

where f = δTi/T0 or Viy, t1 = 140 s and t2 = 500 s. For all studied periods, the maximum heating rate was observed below y = 1.2 Mm, which agrees with Fig. 3. Additionally, with decreasing driver period, ⟨δTi/T0⟩t peaks reach smaller values. Above the altitude y = 1.2 Mm, averaged heating values oscillate at a fixed level. Studies of ⟨Viy⟩t show plasma outflows in the chromosphere.

For comparison purposes, we ran simulations without the SE, corresponding to the case of pure HE. A comparison of Viy (Fig. 6a) with its counterpart for the SE (Fig. 3c) revealed slightly lower plasma velocities in the HE+SE case. This small difference is seen in Fig. 7, which displays Viy for HE (solid line) and HE+SE (dashed line) at the height y = 1.7 Mm and at t = 500 s. The opposite scenario was observed for δTi/T0, which reveals larger values and, thus, more significant chromosphere heating for HE than for HE+SE.

thumbnail Fig. 6.

Time–distance plots for Viy(x = 0, y, t) (a) and δTi/T0 (x = 0, y, t) (b) for the driver parameters corresponding to the panels (c, d) of Fig. 3, but drawn here for HE.

thumbnail Fig. 7.

Variation of Viy(x, y = 1.75 Mm, t = 500 s) s for HE (solid line) and HE+SE (dashed line).

To investigate the influence of ion–neutral collisions in the chromosphere heating, we performed studies of the velocity drift (Fig. 8), Viy − Viy for HE (a) and HE+SE (b). The highest values of Viy − Viy in both cases were observed in the region where the driver was launched, viz. 1 Mm <  y <  1.2 Mm. However, the HE case shows a significant difference in the velocities at the whole chromosphere. These results agree with the heating presented in Figs. 3 and 6. Thus, the heating of the top chromosphere is the result of ion–neutral collisions, which were included in the two-fluid model (Eqs. (1)–(11)). For the HE+SE case, the lowest values of Viy − Vny and δTi/T0 were observed in the middle chromosphere. On the other hand, velocity drift reaches its highest values in the upper chromosphere, which is particularly evident in the HE+SE panel. This is unsurprising as the coupling between species decreases with lower density at higher altitudes. We note that in both equilibrium configurations, the velocities of neutrals are greater than those of ions in most of the chromosphere.

thumbnail Fig. 8.

Time–distance plots for Viy − Vny for the driver parameters corresponding to the panel (c) of Fig. 3, in the case of HE (a) and HE+SE (b).

Figure 9 presents the relative perturbed ion temperature averaged over time and height. This quantity is defined as

H = 1 y 2 − y 1 ∫ y 1 y 2 〈 δ T i T 0 〉 t d y , $$ \begin{aligned} H = \frac{1}{{ y}_{\rm 2} - { y}_{\rm 1}} \int ^{{ y}_{2}}_{{ y}_{\rm 1}} \Big \langle \frac{\delta T_{\rm i}}{T_{\rm 0}} \Big \rangle _t \, \mathrm{d}{ y}, \end{aligned} $$(21)

thumbnail Fig. 9.

Relative perturbed temperature of ions averaged over time and height, H, vs. Pd.

where y1 = 1 Mm and y2 = 1.9 Mm. As H experiences fall-off of with Pd, the chromosphere heating decreases with the driver period. These results converge with that of Soler (2024) (his Eq. (7.1)), which shows that the heating rate is proportional to ω2. Thus, the highest efficiency in heating is achieved with high-frequency waves. A similar effect was observed in the context of Alfvén waves (Song & Vasyliūnas 2011). Thus, we find that the chromosphere acts as a filter, damping high-frequency waves which deposit heat at the lower altitudes.

4. Summary and conclusion

We constructed a new model of the solar chromosphere. In this model, the unperturbed state was described by the magnetohydrostatic equilibrium, with magnetic field lines fanning out with height, and supplemented by the SE to model the variable ionisation degree in the photosphere and chromosphere. This magnetohydrostatic equilibrium is perturbed by a monochromatic driver that operates in the middle chromosphere, exciting essentially SMAWs. This atmosphere is described by two-fluid equations (with ionisation and recombination terms included) for ions (protons) combined with electrons and neutrals (hydrogen atoms), which are solved numerically using the JOANNA code (Wójcik et al. 2019).

We find that, as a result of ion–neutral collisions, the SMAWs are damped and they thermalise their energy at y = 1.05 Mm, with the stronger effect corresponding to smaller wave periods. Most of the energy from granulation-excited waves is concentrated in five-minute oscillations with wave periods of 300 s (Zaqarashvili et al. 2011). Consequently, the 2.5, 5, and 10 s wave periods correspond to less energetic oscillations. These shorter periods are closer to the timescales of ion–neutral collisions, making the two-fluid effects more significant for these waves. Velocity drift reaches the highest values at the bottom of the numerical region, showing agreement with the heating profile. Additionally, the SMAW drive flows from the middle chromosphere into the solar corona. Thus, we find that the realistic ionisation rate specified using the SE results in lower plasma outflows and chromosphere heating than the ionisation determined by the pure HE. We therefore conclude that the SE plays an important role in the evolution of SMAWs. Comparison with observational radiative energy loss data shows that the kinetic energy flux values we obtained are too small to drive the solar wind.

Acknowledgments

KM’s work was done within the framework of the project from the Polish Science Center (NCN) Grant No. 2020/37/B/ST9/00184. Dr. Luis Kadowaki implemented the HE+SE equilibrium into the JOANNA code, funded by this grant. A part of numerical simulations was run on the LUNAR cluster at the Institute of Mathematics at M. Curie-Skłodowska University in Lublin, Poland. The authors gratefully acknowledge Poland’s high-performance computing infrastructure, PLGrid (HPC Centers: ACK Cyfronet AGH), for providing computer facilities and support within computational grant no. PLG/2022/015868. The simulation data was visualised using Python scripts. SP acknowledges support from the projects C14/19/089 (C1 project Internal Funds KU Leuven), G0B5823N and G002523N (WEAVE) (FWO-Vlaanderen), 4000145223 (SIDC Data Exploitation (SIDEX2), ESA Prodex), and Belspo project B2/191/P1/SWiM, as well as from SWATNet, a project that has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 955620.

References

  1. Avrett, E. H., & Loeser, R. 2008, ApJS, 175, 229 [Google Scholar]
  2. Banerjee, D., Krishna Prasad, S., Pant, V., et al. 2021, Space Sci. Rev., 217, 76 [NASA ADS] [CrossRef] [Google Scholar]
  3. Biermann, L. 1946, Naturwissenschaften, 33, 118 [Google Scholar]
  4. Biermann, L. 1947, Naturwissenschaften, 34, 87 [NASA ADS] [CrossRef] [Google Scholar]
  5. Botha, G. J. J., Arber, T. D., Nakariakov, V. M., & Zhugzhda, Y. D. 2011, ApJ, 728, 84 [Google Scholar]
  6. Braginskii, S. I. 1965, Rev. Plasma Phys., 1, 205 [Google Scholar]
  7. Defouw, R. J. 1976, ApJ, 209, 266 [NASA ADS] [CrossRef] [Google Scholar]
  8. Duckenfield, T. J., Kolotkov, D. Y., & Nakariakov, V. M. 2021, A&A, 646, A155 [EDP Sciences] [Google Scholar]
  9. Durran, D. R. 2010, Numerical Methods for Fluid Dynamics: With Applications to Geophysics (Springer Science& Business Media), 32 [CrossRef] [Google Scholar]
  10. Felipe, T., Khomenko, E., & Collados, M. 2010, ArXiv e-prints [arXiv:1005.3684] [Google Scholar]
  11. Felipe, T., Kuckein, C., & Thaler, I. 2018, A&A, 617, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. González-Morales, P. A., Khomenko, E., Vitas, N., & Collados, M. 2020, A&A, 642, A220 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Khomenko, E., & Collados, M. 2015, Liv. Rev. Sol. Phys., 12, 6 [Google Scholar]
  14. Khomenko, E., Collados, M., Díaz, A., & Vitas, N. 2014, Phys. Plasmas, 21, 092901 [Google Scholar]
  15. Kraśkiewicz, J., Murawski, K., & Musielak, Z. E. 2019, A&A, 623, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Kraśkiewicz, J., Murawski, K., & Musielak, Z.E. 2022, MNRAS [arXiv:2211.16463] [Google Scholar]
  17. Kraśkiewicz, J., Murawski, K., Zhang, F., & Poedts, S. 2023, Sol. Phys., 298, 11 [CrossRef] [Google Scholar]
  18. Krishna Prasad, S., Jess, D. B., & Khomenko, E. 2015, ApJ, 812, L15 [NASA ADS] [CrossRef] [Google Scholar]
  19. Kuźma, B., & Murawski, K. 2018, ApJ, 866, 50 [CrossRef] [Google Scholar]
  20. Kuźma, B., Murawski, K., Musielak, Z. E., Poedts, S., & Wójcik, D. 2021, A&A, 652, A88 [EDP Sciences] [Google Scholar]
  21. MacBride, C. D., Jess, D. B., Khomenko, E., & Grant, S. D. T. 2022, ApJ, 938, 154 [NASA ADS] [CrossRef] [Google Scholar]
  22. Maneva, Y. G., Alvarez Laguna, A., Lani, A., & Poedts, S. 2017, ApJ, 836, 197 [Google Scholar]
  23. Miyamoto, S. 1949, PASJ, 1, 10 [NASA ADS] [Google Scholar]
  24. Miyoshi, T., & Kusano, K. 2005, J. Comput. Phys., 208, 315 [NASA ADS] [CrossRef] [Google Scholar]
  25. Murawski, K., Musielak, Z. E., & Wójcik, D. 2020, ApJ, 896, L1 [NASA ADS] [CrossRef] [Google Scholar]
  26. Murawski, K., Musielak, Z. E., Poedts, S., Srivastava, A. K., & Kadowaki, L. 2022, Ap&SS, 367, 111 [NASA ADS] [CrossRef] [Google Scholar]
  27. Musielak, Z. E., Musielak, D. E., & Mobashi, H. 2006, Phys. Rev. E, 73, 036612 [NASA ADS] [CrossRef] [Google Scholar]
  28. Nakariakov, V. M., Afanasyev, A. N., Kumar, S., & Moon, Y. J. 2017, ApJ, 849, 62 [Google Scholar]
  29. Nakariakov, V. M., Kosak, M. K., Kolotkov, D. Y., et al. 2019, ApJ, 874, L1 [Google Scholar]
  30. Niedziela, R., Murawski, K., & Poedts, S. 2021, A&A, 652, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Niedziela, R., Murawski, K., Kadowaki, L., Zaqarashvili, T., & Poedts, S. 2022, A&A, 668, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Ofman, L. 2010, Liv. Rev. Sol. Phys., 7, 4 [Google Scholar]
  33. Ofman, L., & Wang, T. 2022, ApJ, 926, 64 [NASA ADS] [CrossRef] [Google Scholar]
  34. Oliver, R., Soler, R., Terradas, J., & Zaqarashvili, T. V. 2016, ApJ, 818, 128 [Google Scholar]
  35. Pelekhata, M., Murawski, K., & Poedts, S. 2021, A&A, 652, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Pelekhata, M., Murawski, K., & Poedts, S. 2022, ArXiv e-prints [arXiv:2211.12898] [Google Scholar]
  37. Popescu Braileanu, B., Lukin, V. S., Khomenko, E., & de Vicente, Á. 2019, A&A, 627, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Saha, M. N. 1920, Nature, 105, 232 [CrossRef] [Google Scholar]
  39. Schwarzschild, M. 1948, ApJ, 107, 1 [Google Scholar]
  40. Soler, R. 2024, Phil. Trans. R. Soc. London Ser. A, 382, 20230223 [Google Scholar]
  41. Song, P., & Vasyliūnas, V. M. 2011, J. Geophys. Res. (Space Phys.), 116, A09104 [NASA ADS] [Google Scholar]
  42. Stark, B. A., & Musielak, Z. E. 1993, ApJ, 409, 450 [Google Scholar]
  43. Suzuki, T. K., Shimizu, K., & Shoda, M. 2022, Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun (), 147 [Google Scholar]
  44. Vranjes, J., & Krstic, P. S. 2013, A&A, 554, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Withbroe, G. L., & Noyes, R. W. 1977, ARA&A, 15, 363 [Google Scholar]
  46. Wójcik, D., Kuźma, B., Murawski, K., & Srivastava, A. K. 2019, ApJ, 884, 127 [CrossRef] [Google Scholar]
  47. Wójcik, D., Kuźma, B., Murawski, K., & Musielak, Z. E. 2020, A&A, 635, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Zaqarashvili, T. V., Murawski, K., Khodachenko, M. L., & Lee, D. 2011, A&A, 529, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

All Figures

thumbnail Fig. 1.

Variation with height of the equilibrium temperature, T0 (a), ionisation degree, I (b) for HE (dashed line) and HE+SE (solid line), and mass density (c, d) for ions (solid line) and neutrals (dashed line) for HE+SE (c) and HE (d).
In the text
thumbnail Fig. 2.

Spatial profile of Vi(x, y) (expressed in km s−1) at t = 110 s overlaid by its vectors (white arrows) and magnetic field lines for Pd = 5 s.
In the text
thumbnail Fig. 3.

Time-distance plots for Viy(x = 0, y, t) (a, c, e) and δTi/T0 (x = 0, y, t) (b, d, f) for Pd = 10 s (a, b), Pd = 5 s (c, d), and Pd = 2.5 s (e, f), in the case of HE+SE.
In the text
thumbnail Fig. 4.

Time–distance plots (a, b) and averaged-over-time vertical profiles (c, d) for the kinetic energy flux FEics (a, c) and frictional heating term Q (b, d) for Pd = 5  s, in the case of HE+SE.
In the text
thumbnail Fig. 5.

Temporarily averaged ⟨δTi/T0⟩t (a, c, e) and ⟨Viy⟩t (b, d, f) for Pd = 10 s (a, b), Pd = 5 s (c, d), and Pd = 2.5 s (e, f), in the case of HE+SE.
In the text
thumbnail Fig. 6.

Time–distance plots for Viy(x = 0, y, t) (a) and δTi/T0 (x = 0, y, t) (b) for the driver parameters corresponding to the panels (c, d) of Fig. 3, but drawn here for HE.
In the text
thumbnail Fig. 7.

Variation of Viy(x, y = 1.75 Mm, t = 500 s) s for HE (solid line) and HE+SE (dashed line).
In the text
thumbnail Fig. 8.

Time–distance plots for Viy − Vny for the driver parameters corresponding to the panel (c) of Fig. 3, in the case of HE (a) and HE+SE (b).
In the text
thumbnail Fig. 9.

Relative perturbed temperature of ions averaged over time and height, H, vs. Pd.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.