Issue 
A&A
Volume 639, July 2020



Article Number  A45  
Number of page(s)  10  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/201937260  
Published online  06 July 2020 
Numerical simulations of the lower solar atmosphere heating by twofluid nonlinear Alfvén waves
^{1}
Institute of Physics, University of Maria CurieSkłodowska, Pl. M. CurieSkłodowskiej 5, 20031 Lublin, Poland
^{2}
Institute of Space Science and Applied Technology, Harbin Institute of Technology, Shenzhen, Guangdong 518055, PR China
email: yuanding@hit.edu.cn
^{3}
Center for Mathematical Plasma Astrophysics, Department of Mathematics, KU Leuven, Celestijnenlaan 200B, 3001 Leuven, Belgium
Received:
5
December
2019
Accepted:
19
May
2020
Context. We present new insight into the longstanding problem of plasma heating in the lower solar atmosphere in terms of collisional dissipation caused by twofluid Alfvén waves.
Aims. Using numerical simulations, we study Alfvén wave propagation and dissipation in a magnetic flux tube and their heating effect.
Methods. We set up 2.5dimensional numerical simulations with a semiempirical model of a stratified solar atmosphere and a forcefree magnetic field mimicking a magnetic flux tube. We consider a partially ionized plasma consisting of ion + electron and neutral fluids, which are coupled by ionneutral collisions.
Results. We find that Alfvén waves, which are directly generated by a monochromatic driver at the bottom of the photosphere, experience strong damping. Lowamplitude waves do not thermalize sufficient wave energy to heat the solar atmospheric plasma. However, Alfvén waves with amplitudes greater than 0.1 km s^{−1} drive through ponderomotive force magnetoacoustic waves in higher atmospheric layers. These waves are damped by ionneutral collisions, and the thermal energy released in this process leads to heating of the upper photosphere and the chromosphere.
Conclusions. We infer that, as a result of ionneutral collisions, the energy carried initially by Alfvén waves is thermalized in the upper photosphere and the chromosphere, and the corresponding heating rate is large enough to compensate radiative and thermalconduction energy losses therein.
Key words: Sun: activity / Sun: chromosphere / Sun: transition region / sunspots / magnetohydrodynamics (MHD) / waves
© ESO 2020
1. Introduction
A key unresolved issue in solar physics is the energy transport and dissipation from the lower atmosphere to the higher layers of the Sun. The chromospheric plasma loses its thermal energy more rapidly than the corona, which, according to our current understanding, would not be compensated sufficiently by energy coming from lower regions. Therefore, an additional source of chromospheric heating is required (see e.g., Narain & Ulmschneider 1996). The idea that energy carried by waves can heat the solar atmosphere was developed by Bierman (1946) and Schwarzschild (1948), who originally proposed that acoustic waves could be the main heating agents of the chromosphere, with their wavefronts steepening with height, and energy dissipating in the shocks grown in the chromosphere. The photospheric plasma exhibits a wide range of fluid dynamics and is considered to be the main source of waves in the upper atmosphere. Among diverse waves, Alfvén waves can be excited there and propagate into the upper solar atmosphere.
Alfvén waves were theoretically predicted by Alfvén (1942) and subsequently their existence in laboratory plasma was confirmed by Lundquist (1949). These waves were also observed in natural plasmas such as aurora (Chmyrev et al. 1988), the solar wind (Belcher & Davis 1971), and potentially in the solar atmosphere (Tomczyk et al. 2007), although the latter was debated as the observed waves could be interpreted as fast magnetoacoustic kink waves (e.g., Doorsselaere et al. 2008). On the other hand, recent observations suggest the appearance of shortperiod torsional Alfvén waves propagating in finely structured flux tubes (Srivastava et al. 2017). These structures were further investigated numerically by Kuźma & Murawski (2018), who showed that finely structured flux tubes act as wave guides for both torsional Alfvén and kink waves and that the latter transport most of the energy flux into the corona.
So far, the wave processes in the solar atmosphere have mainly been investigated within the framework of magnetohydrodynamics (MHD) as a valid model for a fully ionized plasma. However, the relatively low plasma temperatures in the photosphere and chromosphere result in partial ionization of the plasma. Therefore, both these regions of the atmosphere should be modeled while taking into account the presence of neutrals (e.g., Zaqarashvili et al. 2011). Efforts in this direction were made by MartínezSykora et al. (2012) and Khomenko (2017) who solved numerically the 2.5D MHD equations with ambipolar diffusion and showed that Alfvén waves are damped effectively there. The minor ions collide frequently with the neutrals in low atmospheric regions and lead to significant damping of Alfvén waves (Piddington 1956; Kulsrud & Pearce 1969; Vranjes & Poedts 2009, 2010; Soler et al. 2014); this damping and the associated thermal energy deposition can compensate the radiative losses there (Song & Vasyliūnas 2011; Goodman 2011). Zaqarashvili et al. (2013) showed analytically that highfrequency Alfvén waves are damped efficiently by ionneutral collisions, whereas Murawski & Musielak (2010) and Wójcik et al. (2017) proposed that lowfrequency Alfvén waves are evanescent as a result of the solar gravitational stratification. Vranjes & Poedts (2010) suggested that driftAlfvén waves could heat the corona stochastically.
Recent work on twofluid nonlinear Alfvén waves in the solar atmosphere by MartínezGómez et al. (2018) revealed frictional heating due to ionneutral collisions and a consequent rise of the plasma temperature. As an application, these authors investigated physical conditions akin to those in a solar quiescent prominence. They concluded also that the nonlinear effects are more pronounced for standing oscillations than for propagating waves. Soler et al. (2019) discussed heating of partially ionized plasma in an expanding magnetic flux tube embedded in the stratified lower atmosphere. These latter authors considered the steady state of wave propagation, implying a driver that acts continuously and for a sufficiently long time for a stationary Alfvén wave pattern to be achieved. We extend the previously devised models by taking into account a realistic Alfvén wave propagation in a stratified twofluid medium of magnetic flux tube. We examine effects of Alfvén wave reflection, ionneutral collisional damping, and excitation of magnetoacoustic waves via the ponderomotive force.
We organize this paper as follows. In Sect. 2 we present the twofluid equations used in this study and the numerical implementation. In Sect. 3 we present our model of the stratified atmosphere and magnetic field, and in Sect. 4 we show and discuss the results of our numerical simulations. Finally, in Sect. 5 we summarize and present our conclusions.
2. Twofluid plasma model
In this study we take into account hydrogen as the main plasma constituent and model the influence of heavier elements by use of the OPAL repository of solar abundances (e.g., Vögler et al. 2005). Additionally, we assume that ions and electrons constitute a single ion–electron fluid, whereas neutrals are described by another fluid that interacts with the ions via ionneutral collisions. In the model we neglect several other effects such as recombination and ionization, charge exchange and the effects of electrons exerted on ions, which can be modeled with extra terms in the generalized Ohm’s law (e.g., Khomenko 2017).
We describe the solar atmosphere by a twofluid model including the mass continuity, momentum, and energy equations. The ion and electron version of these equations is supplemented by the induction equation in which the effect of electrons is neglected because their mass is significantly smaller than the mass of ions and neutrals (e.g., Ballester et al. 2018a and references therein). The hydrodynamic (HD) equations for the neutrals are:
and the MHD equations for the ions + electrons are:
denote the total energy densities of neutrals and ions, respectively. These equations are supplemented by the ideal gas laws for both neutrals and ions:
Here, 𝜚_{n} (𝜚_{i}) is the mass density of the neutrals (ions), V_{n} (V_{i}) the velocity of the neutrals (ions), p_{n} (p_{ie}) the neutral (ion and electron) gas pressure, B the magnetic field, μ the magnetic permeability of the medium, γ = 5/3 the adiabatic index, and g = [0, −g, 0] is the vector of gravitational acceleration with its magnitude g = 274.78 m s^{−2} acting along yaxis in Cartesian coordinates. In addition, T_{i} (T_{n}) is the ion (neutral) temperature, μ_{ie} = 0.58 and μ_{n} = 1.21 are the mean masses of respectively electrons + ions and neutrals, m_{H} is the hydrogen mass, k_{B} is Boltzmann’s constant, and α_{in} is ionneutral drag coefficient linked to ionneutral collision crosssection, σ_{in}, by the following formula (e.g., Oliver et al. 2016; Ballester et al. 2018a):
Here we take a quantum value of σ_{in} (Vranjes & Krstic 2013). Additionally, Q_{i} and Q_{n} denote the heat production and exchange terms resulting from ionneutral collisions such as (Ballester et al. 2018a and references therein):
This additional heating results solely from frictional interaction between both fluids and is proportional to the velocity difference squared between both species and the friction coefficient α_{in}. The second terms in Eqs. (11) and (12) denote heat exchange between ions and neutrals.
Henceforth, we limit our discussion to a 2.5D Cartesian region with x and y corresponding to horizontal and vertical axes, respectively. We consider a region spanning from the photosphere up to the low corona. We assume that the system is invariant along the zdirection with ∂/∂z = 0 but the perturbations of the velocities (V_{i z} and V_{n z}) and magnetic field (B_{z}) components along this direction are not identically zero. In such geometry, V_{i z} and B_{z} correspond to Alfvén waves which decouple from the magnetoacoustic waves in this setup, while V_{n z} mimics the strongly damped neutral entropy waves (Zaqarashvili et al. 2011, Soler et al. 2013).
3. Hydrostatic equilibrium model of the solar atmosphere
In this section, we consider the equilibrium conditions of the solar atmosphere. We assume that this atmosphere is initially in hydrostatic equilibrium.
3.1. Ion and neutral mass densities
The hydrostatic equilibrium gas pressure and mass density (Fig. 1, top) profiles are specified uniquely by the temperature profile, which we take here to be the same for ions and neutrals (e.g., Oliver et al. 2016) from the quiet solar atmosphere semiempirical model of Avrett & Loeser (2008) (Fig. 1, bottom). Using ideal gas law (Eq. (9)) and hydrostatic (−∇p_{ie, n} + 𝜚_{i, n} g = 0) equations we obtain:
Fig. 1. Vertical profiles of the equilibrium ion, 𝜚_{i}(y), and neutral mass densities, 𝜚_{n}(y), (top; dotted and dashed lines respectively), and the hydrostatic equilibrium temperature, T(y), taken from Avrett & Loeser (2008; bottom). 
where p_{0 n} (p_{0 i}) is neutral (ion) pressure at the reference point y_{ref} = 50 Mm, and
are respectively the neutral and ion pressure scaleheights.
As a consequence of the temperature variation, the ionization degree, specified as 𝜚_{i}(y)/𝜚_{n}(y), depends on altitude in the solar atmosphere. In the photosphere, which is the layer occupying the region 0 < y < 0.5 Mm and where the temperature is only about 5600 K, the plasma is weakly ionized with only one ion for every approximately 10^{3} neutrals, while in the hot solar corona (∼1 million K), which spans from the transition region located at y ≈ 2.1 Mm (see Fig. 1, bottom), the plasma is essentially fully ionized. Figure 1 (top) illustrates the vertical profiles of the mass densities of ions and neutrals. We note that the ion mass density is much lower than the neutral mass density in the photosphere and the low chromosphere. These two densities become comparable about 900km below the transition region. Higher up, the neutrals are more abundant than ions and in the low corona the mass density of the neutrals experiences a sudden falloff with height.
In our hydrostatic equilibrium, the temperature of both ions and neutrals (Fig. 1, bottom) reaches a minimum of T = 4300 K at the top of the photosphere (500km) and increases suddenly in the transition region between the chromosphere and the low corona, which is located at y = 2.1 Mm.
3.2. Currentfree magnetic field
In the equilibrium model, we adopt a currentfree (∇ × B = 0) magnetic field, B = [B_{x}, B_{y}, B_{z}], which describes a solar magnetic fluxtube (Low 1985):
where B_{v} = 5 G is the magnitude of a straight vertical magnetic field component, a = −0.3 Mm is the vertical location of the singularity in the magnetic field, and S is a free parameter corresponding to the magnetic field strength. We set S so that the maximum value of B ≈ 82 G and is reached at the fluxtube foot point, (x = 0, y = 0) Mm. The spatial profile of the magnetic field lines and the corresponding magnetic field strength, B(x, y), are displayed in Fig. 2 (left). We note that magnitude of the magnetic field falls off from the fluxtube foot point and the magnetic field lines become essentially vertical in the corona.
Fig. 2. Spatial profiles of equilibrium magnetic field lines (left) and the Alfvén speed, c_{A}(x, y), (right). The magnetic field strength, B(x, y), is expressed in Gauss, while c_{A} is expressed in km s^{−1}. Only a part of the right side of the simulation box is displayed. 
The spatial profile of the Alfvén speed, , is illustrated in Fig. 2 (right). In the upper photosphere and the lower chromosphere, the Alfvén speed c_{A} attains values lower than 10 km s^{−1} along the fluxtube center, given by x = 0 up to y ≈ 0.6 Mm. Higher up, as a result of the low total ion + neutral mass density, c_{A} rises to values above 40 km s^{−1} until the transition region, where it experiences a sudden growth to values close to 300 km s^{−1}.
3.3. Analytically estimated wave damping
In this section we first discuss the results of analytical estimates of the damping rate of Alfvén waves in the photosphere and chromosphere.
The dispersion relation for Alfvén waves propagating in a homogeneous medium can be written as (Zaqarashvili et al. 2011; Soler et al. 2013; Ballester et al. 2018a)
where
are respectively the ionneutral and neutralion collision frequencies, χ = 𝜚_{n}/𝜚_{i} denotes ionization fraction, k is the wave number, ω the cyclic frequency, and i stands for the imaginary unit, such that i^{2} = −1. We solve the above dispersion relation for a real value of ω = 2π/P_{d}, where P_{d} is the wave period, to find the real, k_{r}, and imaginary, k_{i}, parts of k, given as
The above result allows us to estimate the wavelength, λ_{r} = 2π/k_{r}, and the damping length, λ_{i} = 1/k_{i}, of the excited Alfvén waves.
Figure 3 illustrates the analytically obtained values of λ_{i} for both photospheric (solid line) and chromospheric (dashed line) conditions. The photosphere is dominated by neutrals, and therefore the specified values of ion and neutral mass densities are respectively 𝜚_{i} = 4.7 × 10^{−10} g cm^{−3} and 𝜚_{n} = 2 × 10^{−7} g cm^{−3}. The corresponding values in the chromosphere are 𝜚_{i} = 1.7 × 10^{−12} g cm^{−3} and 𝜚_{n} = 2.1 × 10^{−12} g cm^{−3}, and as a result χ in the photosphere and chromosphere differs by three orders of magnitude. A strongly coupled photospheric plasma is characterized by a friction coefficient α_{in} seven orders of magnitude larger there than in the chromosphere (6.5 × 10^{−3} vs. 2.6 × 10^{−10}). The magnetic flux tube, given by Eq. (15), results in c_{A} = 3.1 km s^{−1} at the photospheric level and c_{A} = 30 km s^{−1} in the chromosphere. We conclude that Alfvén waves of P_{d} < 0.1 s are strongly damped in the photosphere on a distance spanning from about a few meters up to 15 km. The situation is different in the chromosphere, however. Although wave periods shorter than P_{d} ≈ 0.1 s are still strongly damped there, the damping length is almost constant at a value of λ_{i} ≈ 2.3 km. Waves with longer wave periods are more weakly damped and propagate up to the transition region and higher up into the corona. We note that for P_{d} > 1, the damping lengths in the photosphere and chromosphere are identical. This scenario may differ depending on the model. Cranmer & van Ballegooijen (2005) investigated the process of Alfvén wave reflection in the stratified atmosphere, concluding that Alfvén waves are reflected mostly below the transition region and this process is at its maximum at the photosphere with the reflection coefficient of 0.99825 for a period of about 40 min, whereas shorter wave periods experience weaker reflection. In our simulations, Alfvén waves are affected by reflections in the stratified atmosphere, essentially from regions of abrupt variations in Alfvén speed. However, the wave damping remains the dominant process in the attenuation of Alfvén wave amplitude.
Fig. 3. Analytically obtained damping length, λ_{i}, of Alfvén waves in the photosphere (solid line) and chromosphere (dashed line) vs. wave period P_{d}. 
4. Numerical results
We use the JOANNA code to solve Eqs. (1)–(7). This code is based on highorder Godunovtype methods (Wójcik et al. 2018, 2019) and uses the divergencecleaning method of Dedner et al. (2002) to control the ∇ ⋅ B = 0 constraint.
Unless otherwise stated, we extend the 2.5D simulation box from −2.56 Mm to 2.56 Mm in the xdirection and from 0 to 50 Mm in the ydirection, and cover the simulation region above the photosphere (y = 0) and below the low corona (y = 2.56 Mm) by a uniform grid of 1024 × 512 cells which leads to a uniform grid with a spatial resolution of Δx = Δy = 5 km. We chose this grid to minimize numerical diffusion which essentially does not affect Alfvén wave propagation from the lower atmospheric regions. Above the level y = 2.56 Mm we divide the simulation box into several patches of progressively larger cells along the ydirection. Such a stretched grid plays the role of a sponge as it absorbs the incoming signal and minimizes spurious reflections from the top boundary which could pollute our simulations.
We implement fixed boundary conditions at the bottom and top of the simulation box at which we set all plasma quantities to their equilibrium values. The only exception are the transverse components of the ion and neutral velocities, which in order to study Alfvén waves are specified in the form of a periodic driver acting at the bottom of the photosphere at y = 0. Specifically, we define
where V_{0} is the amplitude of the driver, w its width, and P_{d} the period. We allow V_{0} to vary in the way that the effective amplitude of the driver lies within the range of 0.01 km s^{−1} and 1 km s^{−1}. These amplitudes are in the observational context of a powerlaw spectrum for the quiet Sun. EUV intensity and Doppler shift variations (Mathioudakis et al. 2013) for lowenergy and shortwave periods, could mimic the action of pmodes. However, the latter exhibit much lower amplitudes (≈10^{−4} km s^{−1}) than those used in our work. Instead, our periodic driver is associated with the solar granulation with horizontal velocities of a fraction of 1 km s^{−1}. Additionally, we impose open boundary conditions for outflowing signals at the side boundaries, x = ±2.56 Mm.
4.1. Numerically estimated Alfvén waves damping
In this section we discuss the numerically obtained results about the propagation and damping rate of Alfvén waves in this twofluid setup.
First, we run the 1D simulations of damping of Alfvén waves in the photosphere. We set our numerical experiment with a uniform, highspatialresolution grid of Δx = Δy = 50 m, and uniform vertical magnetic field with its magnitude B = 50 Gauss. Such a fine grid allows us to significantly reduce numerical damping of the waves in our system. We specify the driver at the bottom boundary by Eq. (19), and run these simulations for a wide range of amplitudes, starting with the linear case of 1 m s^{−1} up to a strongly nonlinear case with an amplitude of 40 km s^{−1}. The obtained results reveal strong wave damping that occurs already in the photosphere. In particular, Fig. 4 illustrates the vertical profiles of the zcomponent of the ion velocity, V_{i z}(y), at t = 50 s, t = 200 s, and t = 500 s in the case of the driving period P_{d} = 30 s and the amplitude V_{0} = 40 km s^{−1}, that is, the strongly nonlinear case. The Alfvén waves are completely attenuated before they can reach the height y = 0.15 Mm. We note that this attenuation length does not correspond to the collisional damping length resulting from Eq. (16). For P_{d} < 1 s, wave damping is due to ion–neutral collisions, while for P_{d} > 1 s, the Alfvén waves are affected by the cutoff period and the turning point, the latter being very low for longperiod waves (Wójcik et al. 2017). Further investigation of nonlinear Alfvén waves will be the subject of future studies.
Fig. 4. Vertical profiles of the zcomponent of the ion velocity, V_{i z}(y), at t = 50 s (solid line), t = 200 s (dashed line), and t = 500 s (dotted line) in the case of the driving period P_{d} = 30 s and the amplitude V_{0} = 40 km s^{−1}, for the uniform vertical magnetic field of B_{y} = 50 G. 
Figure 5 shows the spatial profile of the zcomponent of the ion velocity, V_{i z}(x, y, t), obtained by means of numerical simulations for the magnetic field given by Eq. (15) at six instants in time, namely at t = 180 s, t = 250 s, t = 500 s, t = 750 s, t = 1250 s, and t = 2300 s (from topleft to bottomright) for P_{d} = 30 s and V_{0} = 10 m s^{−1}. We clearly see that the driver excites transversal plasma oscillations propagating along the expanding solar magnetic flux tube. These waves are strongly damped in the upper photosphere (0 < y < 0.3 Mm). Higher up, in the chromosphere, the amplitude exponential growth with the distance equal to the pressurescale height (e.g., Kuźma & Murawski 2018) becomes dominant. At these altitudes the coupling between ions and neutrals is already weaker, and therefore neutrals do not follow exactly the same evolution as the ions (not shown).
Fig. 5. Spatial profiles of the zcomponent of the ion velocity, V_{i z}(x, y), at t = 180 s, t = 250 s, t = 500 s, t = 750 s, t = 1250 s, and t = 2300 s (from topleft to bottomright) expressed in units of 1 m s^{−1}. The case of the driving period P_{d} = 30 s and the amplitude V_{0} = 10 m s^{−1}. 
The abovementioned damping of Alfvén waves results from the fact that the transversal neutral velocity component, V_{n z}, is strongly damped (e.g., Zaqarashvili et al. 2011). This results in a significant difference between the transversal components of the ion and neutral velocities, V_{i z} − V_{n z}. Figure 6 (topleft) illustrates the heating rate resulting from ionneutral collisions, Q_{i}, which following Eq. (11) coincides with the ion–neutral drift, V_{i} − V_{n} (not shown). We infer that due to a large value of this drift, a fraction of the energy initially carried by Alfvén waves is thermalized in the bottom layer of the photosphere. Indeed, Fig. 6 (topright) reveals a growth of the plasma temperature close to the bottom boundary up to δT_{i} ≈ 200 K over a timescale of 250 s. We note that acoustic (Kuźma et al. 2019) and fast magnetoacoustic waves (Popescu Braileanu et al. 2019) behave differently from Alfvén waves, as they contribute to plasma heating in the chromosphere, whereas Alfvén waves directly thermalize their energy in the low photosphere. A first attempt to understand how the heating and cooling processes modify the properties of MHD waves in a fully ionized medium in the case of solar prominence was made by Ballester et al. (2016). These studies were later extended to the partially ionized case (Ballester et al. 2018b), and the authors revealed that wave periods and damping times of MHD waves become time dependent when considering the heating and cooling processes. Moreover, in the case of Alfvén waves, the cutoff wave periods also become time dependent and the attenuation rate is completely different. Moreover, this temperature gradient in the solar atmosphere would be subject to thermal conduction, which is not considered in our model. Indeed, Soler et al. (2019) suggested, that thermal conduction and radiative loses could play an influential role in Alfvénwave propagation and damping in the solar atmosphere and that such simulations should be performed. We are considering this as a topic for a future study.
Fig. 6. Timedistance plots of the heating rate resulting from ion–neutral collisions, Q_{i}(y, t) (topleft), perturbed ion temperature, δT_{i}(y, t) (topright), the squared transversal component of the magnetic field, (bottomleft), and the vertical component of ion velocity, V_{i y} (bottomright) evaluated at x = 0.1 Mm for the driving period P_{d} = 30 s and the amplitude V_{0} = 0.4 km s^{−1}. 
However, Alfvén waves drive magnetoacoustic waves through a gradient in magnetic pressure. This is illustrated in Fig. 6 (bottomleft) and means that Alfvén waves act as a driver generating magnetoacoustic waves. These waves are associated with V_{i y} and they propagate upwards into higher layers of the chromosphere (bottomright). Hence, Alfvén waves, which are damped in the photosphere, lead to both a local temperature increase at the bottom of the photosphere and the excitation of magnetoacoustic waves.
The amplitude of the magnetoacoustic waves that are driven by Alfvén waves) is proportional to the amplitude of the driver. Figure 7 illustrates spatial profiles of for the driving period P_{d} = 30 s and the amplitude V_{0} = 0.4 km s^{−1} at four instants in time, namely at t = 1000 s, t = 1250 s, t = 1500 s, and t = 2000 s (from topleft to bottomright). From this figure, we infer that Alfvénwavedriven magnetoacoustic waves have higher amplitudes than the amplitude of the driver exceeding 1 km s^{−1}. We note that, as in our simulations where we implemented the realistic mass density profile, waves can experience amplitude growth with height without providing any external energy; this growth results from momentum conservation, and the carried energy does not increase with y. We also note that the regions with the largest values of are located in the photosphere (0.1 Mm < y < 0.3 Mm) and in the middle to upper chromosphere (y > 1.25 Mm). In the photosphere, the plasma experiences mostly a drift in horizontal direction, while chromospheric plasma moves essentially along the vertical magnetic field lines. Below, we show that these are the regions where the plasma is heated the most.
Fig. 7. Spatial profiles of , corresponding to magnetoacoustic waves at t = 1000 s, t = 1250 s, t = 1500 s, and t = 2000 s (from topleft to bottomright) expressed in units of 1 km s^{−1}. Here, the driving period P_{d} = 30 s and the amplitude V_{0} = 0.4 km s^{−1}. 
Figure 8 (top row) shows time–distance plots of ion magnetoacoustic waves characterized by perturbations in V_{i x} and V_{i y}. We see that V_{i x} dominates in the upper photosphere (topleft), where the magnetic field lines are strongly curved and almost become horizontal, while V_{i y} prevails in chromospheric regions that are permeated by vertical magnetic field (topright). The energy of these waves is deposited in the form of thermal energy heating the plasma, and the energy is released during ion–neutral collisions; see also Kuźma et al. (2019) who showed that twofluid acoustic waves with their wave periods of 30 − 120 s effectively heat chromospheric plasma. Indeed, Fig. 8 (bottomleft) displays time–distance plots of . The phaseshift between ion and neutral waves in the upper photosphere and the chromosphere results in a significant difference between and . From Eq. (11) we infer that this velocity difference results in nonzero values of Q_{i} and thus in a growth of δT_{i} in time. Figure 8 (bottomright) displays a time–distance plot of the perturbed ion temperature averaged over the whole horizontal distance (d = 5.12 Mm), . Analyzing ⟨δT_{i}(y)⟩ we find that plasma in the lower photosphere 0 < y < 0.3 Mm and in the chromosphere 1 Mm < y < 1.5 Mm is significantly heated. This heating coincides with V_{i x} in the photosphere and V_{i y} in the chromosphere. We note that the excited magnetoacoustic waves have wave periods close to 180 s. This periodicity may result from local a cutoff period and all generated wave periods will tend to this value; see also Wójcik et al. (2017) for the discussion on Alfvén wave cutoffs.
Fig. 8. Time–distance plots of V_{i x}(x = 0.1 Mm,y, t) (topleft) and V_{i y}(x = 0.1 Mm,y, t) (topright), Mm, y, t) (bottomleft) expressed in units of 1 km s^{−1}, and ⟨δT_{i}(y, t)⟩ expressed in Kelvin (bottomright) for the driving period P_{d} = 30 s and the amplitude V_{0} = 0.4 km s^{−1}. 
4.2. Heating by twofluid Alfvén waves
In this section we present an estimation of the energy and mass transport by Alfvén waves. The transversal ion mass density and energy fluxes transported by Alfvén waves can be estimated as F_{𝜚 i}(x, y, t) = 𝜚_{i}V_{i z} and (Mathioudakis et al. 2013). The vertical profile of horizontally and temporarily averaged mass and energy fluxes are shown in Fig. 9. Despite the high values of both mass and energy fluxes carried by Alfvén waves in the photosphere, the signal in V_{i z} is attenuated very quickly (solid line). We clearly see that V_{i z} (corresponding to Alfvén waves) is partially transferred into magnetoacoustic waves (denoted by both the dotted line for V_{i x} and by the dashed line representing V_{i y}). Moreover, we conclude that the plasma escaping into the corona above the transition region carries a significant amount of mass density and energy, possibly contributing to overall plasma outflows and solar wind (Wójcik et al. 2019). We note that the obtained values of mass and energy fluxes lay within the range of theoretical findings for energy losses in the upper chromosphere estimated by Withbroe & Noyes (1977). However, this statement does not mean that Alfvén waves effectively contribute to plasma heating. To correctly estimate the impact of propagating waves in terms of plasma heating, we have to look at profiles of ⟨δT_{i}(y)⟩.
Fig. 9. Vertical profiles of horizontally and temporarily averaged ion mass, ⟨F_{𝜚i}(y)⟩, and energy, ⟨F_{E}(y)⟩, fluxes for the driving period P_{d} = 30 s and the amplitude V_{0} = 0.4 km s^{−1}. Solid, dotted, and dashed lines correspond to z, x, and ycomponents, respectively. 
Figure 10 illustrates the heating rate resulting from ion–neutral collisions, that is, the first term in the expression for Q_{i} (see Eq. (11)), versus height, y. This term depends on the ionneutral drift, . The five displayed lines correspond to ion–neutral drift values of 0.01 km s^{−1}, 0.1 km s^{−1}, 1 km s^{−1}, 10 km s^{−1}, and 100 km s^{−1} (from bottom to top), respectively. We infer that in order to effectively heat plasma, the difference between and has to be higher than 1 km s^{−1} or, alternatively, the driver has to work continuously for a very long time. The former requirement is met for steep wavefronts only. Thus, Alfvén waves with amplitudes below 1 km s^{−1} are very unlikely to heat the photosphere and chromosphere, as they do not steepen (enough) with height, while waves with amplitudes V_{0} = 1 km s^{−1} (and larger) can potentially heat the chromosphere.
Fig. 10. Heating rate vs. y. The lines correspond to ion neutral drift, , of 0.01 km s^{−1}, 0.1 km s^{−1}, 1 km s^{−1}, 10 km s^{−1}, and 100 km s^{−1} (from bottom to top). 
Figure 11 shows the maximum of the horizontally and temporarily averaged relative ion perturbed temperature, ⟨δT_{i}⟩(y)/T_{0}(y), versus the driving amplitude V_{0} of Eq. (19) as measured at the top of the photosphere, that is, at y = 0.5 Mm (blue dots) and in the chromosphere at y = 1.25 Mm (red squares). Here T_{0}(y) is the equilibrium temperature (Avrett & Loeser 2008). Intuitively, higher values of V_{0} are supposed to correspond to higher values of ion–neutral drift and as a result lead to higher Q_{i}. Following this pattern, we expect ⟨δT_{i}⟩/T_{0} to grow in the chromosphere with V_{0}. Indeed, Fig. 11 reveals a growth of ⟨δT_{i}⟩/T_{0} with increasing V_{0}. Moreover, we find that Alfvén waves with amplitudes lower than 0.45 km s^{−1} excite magnetoacoustic waves that mostly heat the upper photosphere and the chromosphere. This additional heating is only up to ∼5% of the original photospheric temperature. Alfvén waves of higher amplitudes, namely between 0.5 and 0.7 km ^{−1}, excite magnetoacoustic waves that can penetrate higher layers of the atmosphere and effectively heat the chromosphere up to about twice its initial temperature. Further increase of the driver amplitude, closer to 1 km s^{−1}, leads to plasma heating above ten times its initial temperature. From Fig. 8 (bottomright) we infer that, although the photosphere is heated almost immediately, the chromospheric plasma requires the driver to work continuously for several thousand seconds in order to be heated. This time results from the time needed for Alfvén waves to propagate, to drive magnetoacoustic waves, and for the latter to travel through the chromosphere and reflect from the transition region. We note that after this time the solar granulation in the photosphere is also developed, which can excite magnetoacoustic waves (Wójcik et al. 2019).
Fig. 11. Maximum of horizontally and temporarily averaged relative ion perturbed temperature, ⟨δT_{i}⟩/T_{0}, vs. the driving amplitude V_{0} for the driving period P_{d} = 30 s, collected at y = 0.5 Mm (blue dots) and y = 1.25 Mm (red squares). 
Figure 12 illustrates the dependence of the maximum of the horizontally and temporarily averaged relative ion perturbed temperature, ⟨δT_{i}⟩/T_{0} on the driving period, P_{d}. We infer that the plasma heating rate significantly rises for shortperiod Alfvén waves, while for wave periods longer than P_{d} = 80 s this heating remains essentially constant. The decrease of the heating for longer wave periods is most notably seen in the chromosphere.
Fig. 12. Maximum of horizontally and temporarily averaged relative ion perturbed temperature, ⟨δT_{i}⟩/T_{0}, vs. the driving period P_{d} for the driving amplitude V_{0} = 0.5 km s^{−1}, collected at y = 0.5 Mm (blue dots) and y = 1.25 Mm (red squares). 
5. Conclusions
We devoted this paper to the study of solar atmospheric heating by Alfvén waves in the framework of a twofluid plasma model consisting of ions + electrons and neutrals, where ions are described by the MHD equations and neutrals are governed by the HD equations. These two distinct species of plasma interact with each other by means of collisions modeled by ion–neutral drag forces in the momentum equations.
In our model, Alfvén waves are driven by a monochromatic driver acting at the bottom of the photosphere. We find that due to ion–neutral collisions, Alfvén waves are strongly damped in the photosphere. For small amplitudes of the driver, mainly below V_{0} = 0.1 km s^{−1}, these waves cannot effectively heat the solar atmospheric plasma. Although the obtained values of the energy flux of propagating Alfvén waves lay within the range of the theoretical findings for energy losses in the upper photosphere and chromosphere (Withbroe & Noyes 1977), these lowamplitude waves heat the chromospheric plasma increasing its temperature by only up to 25 K on a timescale of a few minutes. Alfvén waves differ from twofluid acoustic waves which can thermalize their energy in the chromosphere leaving the photospheric plasma essentially intact (Kuźma et al. 2019).
However, Alfvén waves with amplitudes higher than 0.1 km s^{−1}, which can potentially be associated with the solar granulation revealing horizontal velocities being a fraction of 1 km s^{−1}, drive magnetoacoustic waves via the ponderomotive force, and these magnetoacoustic waves can have amplitudes higher than 1 km s^{−1}. Our work confirms the findings of Kuźma et al. (2019) who showed that twofluid acoustic waves heat chromospheric plasma. We find that magnetoacoustic waves driven by nonlinear Alfvén waves can effectively heat both the photosphere and the chromosphere. However, we note that although the photospheric plasma is heated almost immediately, it requires the photospheric driver to act for a sufficiently long time (t > 4000 s) for the chromospheric plasma to be heated in this way.
Further studies are required to determine the possibility of sustaining a quasiequilibrium solar atmosphere heated solely by Alfvén and magnetoacoustic waves, employing 3D simulations and more realistic conditions, including for example radiative losses, thermal conduction, and nonlocal thermal equilibrium in the chromosphere and resistivity terms.
Acknowledgments
This work was done within the framework of the projects from the National Science Centre (NCN) Grant nos. 017/25/B/ST9/00506 and 2017/27/N/ST9/01798. D.Y. was supported by the National Natural Science Foundation of China (NSFC, 11803005, 11911530690), the Shenzhen Technology Project (JCYJ20180306172239618) and Shenzhen Science and Technology program (Group No. KQTD20180410161218820). The JOANNA code was developed by Darek Wójcik. Numerical simulations were performed on the LUNAR cluster at Institute of Mathematics of University of M. CurieSkłodowska, Lublin, Poland.
References
 Alfvén, H. 1942, Nature, 150, 405 [NASA ADS] [CrossRef] [Google Scholar]
 Avrett, E., & Loeser, R. 2008, ApJS, 175, 229 [NASA ADS] [CrossRef] [Google Scholar]
 Ballester, J. L., Carbonell, M., Soler, R., & Terradas, J. 2016, A&A, 591, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ballester, J. L., Alexeev, I., Collados, M., et al. 2018a, Space Sci. Rev., 214, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Ballester, J. L., Carbonell, M., Soler, R., & Terradas, J. 2018b, A&A, 609, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belcher, J. W., & Davis, L. 1971, J. Geophys. Res., 76, 3534 [NASA ADS] [CrossRef] [Google Scholar]
 Bierman, L. 1946, Naturwissenschaften, 33, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Chmyrev, V. M., Bilichenko, S. V., Pokhotelov, O. A., et al. 1988, Phys. Scr., 38, 841 [CrossRef] [Google Scholar]
 Cranmer, S. R., & van Ballegooijen, A. A. 2005, ApJS, 156, 265 [NASA ADS] [CrossRef] [Google Scholar]
 Dedner, A., Kemm, F., Kröner, D., et al. 2002, J. Comput. Phys., 175, 645 [Google Scholar]
 Doorsselaere, T. V., Nakariakov, V. M., & Verwichte, E. 2008, ApJ, 676, L73 [NASA ADS] [CrossRef] [Google Scholar]
 Goodman, M. L. 2011, ApJ, 735, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Khomenko, E. 2017, Plasma Phys. Controlled Fusion, 59, 014038 [NASA ADS] [CrossRef] [Google Scholar]
 Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445 [NASA ADS] [CrossRef] [Google Scholar]
 Kuźma, B., & Murawski, K. 2018, ApJ, 866, 50 [CrossRef] [Google Scholar]
 Kuźma, B., Wójcik, D., & Murawski, K. 2019, ApJ, 878, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Low, B. C. 1985, ApJ, 293, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Lundquist, S. 1949, Phys. Rev., 76, 1805 [CrossRef] [Google Scholar]
 MartínezGómez, D., Soler, R., & Terradas, J. 2018, ApJ, 856, 16 [NASA ADS] [CrossRef] [Google Scholar]
 MartínezSykora, J., Pontieu, B. D., & Hansteen, V. 2012, ApJ, 753, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Mathioudakis, M., Jess, D. B., & Erdélyi, R. 2013, Space Sci. Rev., 175, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Murawski, K., & Musielak, Z. E. 2010, A&A, 518, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Narain, U., & Ulmschneider, P. 1996, Space Sci. Rev., 75, 453 [NASA ADS] [CrossRef] [Google Scholar]
 Oliver, R., Soler, R., Terradas, J., & Zaqarashvili, T. V. 2016, ApJ, 818, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Piddington, J. H. 1956, MNRAS, 116, 314 [NASA ADS] [CrossRef] [Google Scholar]
 Popescu Braileanu, B., Lukin, V. S., Khomenko, E., & de Vicente, Á. 2019, A&A, 627, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schwarzschild, M. 1948, ApJ, 107, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Soler, R., Carbonell, M., & Ballester, J. L. 2013, ApJS, 209, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Soler, R., Ballester, J. L., & Zaqarashvili, T. V. 2014, A&A, 573, A79 [Google Scholar]
 Soler, R., Terradas, J., Oliver, R., & Ballester, J. L. 2019, ApJ, 871, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Song, P., & Vasyliūnas, V. M. 2011, J. Geophys. Res., 116, A09104 [Google Scholar]
 Srivastava, A. K., Shetye, J., Murawski, K., et al. 2017, Sci. Rep., 7, 43147 [NASA ADS] [CrossRef] [Google Scholar]
 Tomczyk, S., McIntosh, S. W., Keil, S. L., et al. 2007, Science, 317, 1192 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vranjes, J., & Krstic, P. S. 2013, A&A, 554, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vranjes, J., & Poedts, S. 2009, MNRAS, 400, 2147 [NASA ADS] [CrossRef] [Google Scholar]
 Vranjes, J., & Poedts, S. 2010, ApJ, 719, 1335 [NASA ADS] [CrossRef] [Google Scholar]
 Withbroe, G. L., & Noyes, R. W. 1977, ARA&A, 15, 363 [NASA ADS] [CrossRef] [Google Scholar]
 Wójcik, D., Murawski, K., Musielak, Z. E., Konkol, P., & Mignone, A. 2017, Sol. Phys., 292, 31 [CrossRef] [Google Scholar]
 Wójcik, D., Murawski, K., & Musielak, Z. E. 2018, MNRAS, 481, 262 [NASA ADS] [CrossRef] [Google Scholar]
 Wójcik, D., Kuźma, B., Murawski, K., & Srivastava, A. K. 2019, ApJ, 884, 127 [CrossRef] [Google Scholar]
 Zaqarashvili, T. V., Khodachenko, M. L., & Rucker, H. O. 2011, A&A, 529, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zaqarashvili, T. V., Khodachenko, M. L., & Soler, R. 2013, A&A, 549, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Figures
Fig. 1. Vertical profiles of the equilibrium ion, 𝜚_{i}(y), and neutral mass densities, 𝜚_{n}(y), (top; dotted and dashed lines respectively), and the hydrostatic equilibrium temperature, T(y), taken from Avrett & Loeser (2008; bottom). 

In the text 
Fig. 2. Spatial profiles of equilibrium magnetic field lines (left) and the Alfvén speed, c_{A}(x, y), (right). The magnetic field strength, B(x, y), is expressed in Gauss, while c_{A} is expressed in km s^{−1}. Only a part of the right side of the simulation box is displayed. 

In the text 
Fig. 3. Analytically obtained damping length, λ_{i}, of Alfvén waves in the photosphere (solid line) and chromosphere (dashed line) vs. wave period P_{d}. 

In the text 
Fig. 4. Vertical profiles of the zcomponent of the ion velocity, V_{i z}(y), at t = 50 s (solid line), t = 200 s (dashed line), and t = 500 s (dotted line) in the case of the driving period P_{d} = 30 s and the amplitude V_{0} = 40 km s^{−1}, for the uniform vertical magnetic field of B_{y} = 50 G. 

In the text 
Fig. 5. Spatial profiles of the zcomponent of the ion velocity, V_{i z}(x, y), at t = 180 s, t = 250 s, t = 500 s, t = 750 s, t = 1250 s, and t = 2300 s (from topleft to bottomright) expressed in units of 1 m s^{−1}. The case of the driving period P_{d} = 30 s and the amplitude V_{0} = 10 m s^{−1}. 

In the text 
Fig. 6. Timedistance plots of the heating rate resulting from ion–neutral collisions, Q_{i}(y, t) (topleft), perturbed ion temperature, δT_{i}(y, t) (topright), the squared transversal component of the magnetic field, (bottomleft), and the vertical component of ion velocity, V_{i y} (bottomright) evaluated at x = 0.1 Mm for the driving period P_{d} = 30 s and the amplitude V_{0} = 0.4 km s^{−1}. 

In the text 
Fig. 7. Spatial profiles of , corresponding to magnetoacoustic waves at t = 1000 s, t = 1250 s, t = 1500 s, and t = 2000 s (from topleft to bottomright) expressed in units of 1 km s^{−1}. Here, the driving period P_{d} = 30 s and the amplitude V_{0} = 0.4 km s^{−1}. 

In the text 
Fig. 8. Time–distance plots of V_{i x}(x = 0.1 Mm,y, t) (topleft) and V_{i y}(x = 0.1 Mm,y, t) (topright), Mm, y, t) (bottomleft) expressed in units of 1 km s^{−1}, and ⟨δT_{i}(y, t)⟩ expressed in Kelvin (bottomright) for the driving period P_{d} = 30 s and the amplitude V_{0} = 0.4 km s^{−1}. 

In the text 
Fig. 9. Vertical profiles of horizontally and temporarily averaged ion mass, ⟨F_{𝜚i}(y)⟩, and energy, ⟨F_{E}(y)⟩, fluxes for the driving period P_{d} = 30 s and the amplitude V_{0} = 0.4 km s^{−1}. Solid, dotted, and dashed lines correspond to z, x, and ycomponents, respectively. 

In the text 
Fig. 10. Heating rate vs. y. The lines correspond to ion neutral drift, , of 0.01 km s^{−1}, 0.1 km s^{−1}, 1 km s^{−1}, 10 km s^{−1}, and 100 km s^{−1} (from bottom to top). 

In the text 
Fig. 11. Maximum of horizontally and temporarily averaged relative ion perturbed temperature, ⟨δT_{i}⟩/T_{0}, vs. the driving amplitude V_{0} for the driving period P_{d} = 30 s, collected at y = 0.5 Mm (blue dots) and y = 1.25 Mm (red squares). 

In the text 
Fig. 12. Maximum of horizontally and temporarily averaged relative ion perturbed temperature, ⟨δT_{i}⟩/T_{0}, vs. the driving period P_{d} for the driving amplitude V_{0} = 0.5 km s^{−1}, collected at y = 0.5 Mm (blue dots) and y = 1.25 Mm (red squares). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.