Monochromatic Two-Fluid Alfvén Waves in the Partially Ionised Solar Chromosphere

We present new results towards the explanation of the chromospheric-heating problem and the solar-wind origin, using a two-fluid model that takes into account the collisional interaction between ions (protons) and neutrals (hydrogen atoms). Our aim is to further reveal the mechanism behind chromospheric heating and plasma outflows. We simulate and analyse the propagation and evolution of Alfvén waves in the partially ionised solar chromosphere, consisting of ions + electrons and neutral fluids. The simplified model chromosphere is permeated by a vertical, uniform magnetic field. We perform numerical simulations in the framework of a quasi-1.5-dimensional (1.5D), two-fluid model in which Alfvén waves are excited by a harmonic driver in the transverse component of the ion and neutral velocities, operating in the chromosphere. In the case of a small-amplitude driver, Alfvén waves are weakly damped, and for the chosen wave periods of a few seconds, Alfvén waves manage to propagate through the chromosphere and enter the solar corona. Non-linear Alfvén waves excited by a large-amplitude driver cause significant chromospheric heating and plasma outflows. We thus conclude that two-fluid Alfvén waves with larger amplitudes can contribute to chromospheric heating and plasma outflows, which may result higher up in the solar-wind origin.


Introduction
The solar atmosphere consists of three distinct layers with different physical properties, and in particular the temperature of the solar atmosphere varies drastically with height. In the bottom layer, called the photosphere, the temperature is only about 5600 K at the solar surface and it decreases to about 4300 K at the top of the photosphere (at about 500 km height), which is not sufficient for a substantial ionisation. In fact, at the top of the photosphere the ionisation degree is as low as 0.01%. Higher up, in the chromosphere (500 -2100 km above the solar surface), the temperature rises to about 6 -7 × 10 3 K at about 1800 km height (and continues to rise faster above this level), and the ionisation rate grows accordingly. In the solar corona, which caps the chromosphere, the temperature rises to about 2 -3 MK on average, and this outer layer of the solar atmosphere is essentially fully ionised. Two of the major outstanding problems in solar physics concern: i) the heating processes that must occur in the solar atmosphere to supply and dissipate enough energy sufficiently rapidly to compensate the energy losses caused by radiation and thermal conduction, and to account for the increase of the temperature in the chromosphere and the narrow transition region between the chromosphere and the corona (Withbroe and Noyes, 1977;Aschwanden, 2005), and ii) the solar-wind origin, i.e. the cause(s) and location of the acceleration of the solarwind plasma (Cranmer, 2019). To understand and explain these two problems, it is worth investigating the energy transport in the lower solar atmosphere, i.e. the photosphere and chromosphere (Wedemeyer-Böhm, Lagg, and Nordlund, 2009).
The photosphere is a dynamical layer and a source of waves propagating into the upper atmosphere. The waves are believed to be responsible for the energy transport in the solar atmosphere. In the present article, we focus on the propagation and dissipation of Alfvén waves. This is a type of transverse magnetohydrodynamic (MHD) wave of ions and electrons that travel along magnetic-field lines with magnetic tension as the restoring force. Such waves have been studied intensively (e.g. Hollweg, 1971;Terradas and Ofman, 2004;Song and Vasyliūnas, 2011;Tu and Song, 2013;Song, 2017) and are observed both in the lower regions of the solar atmosphere (Erdélyi and Fedun, 2007;Jess et al., 2009;Liu et al., 2019;Baker et al., 2021) and in the corona (Tomczyk et al., 2007). Some observations, e.g. Srivastava et al. (2017), indicate that short-period torsional Alfvén waves propagate in finely structured magnetic-flux tubes, which were recently investigated numerically by Kuźma and Murawski (2018). Song (2017) constructed a 2D semi-quantitative model of the chromosphere, suggesting that the Alfvén-wave heating rate and radiative loss are semiquantitatively consistent. In general, Alfvén waves are of great interest as they may participate in the solar atmospheric heating (e.g. Martínez-Gómez, Soler, and Terradas, 2018;Soler et al., 2019) and, thus, in the generation and acceleration of the solar wind (De Pontieu et al., 2007;Hollweg, 2008;Ofman, 2019).
The numerical modelling of Alfvén-wave propagation in the lower solar atmosphere, including the photosphere and chromosphere, is in fact challenging, mainly due to the weakly or partially ionised plasmas involved. Since the single-fluid MHD model is not sufficiently accurate for modelling weakly or partially ionised plasmas that include a significant amount of neutrals, the interactions between ions and neutrals need to be described by ambipolar diffusion (Cally and Khomenko, 2018;Khomenko and Cally, 2019;Nóbrega-Siverio et al., 2020;Popescu Braileanu and Keppens, 2021;MacBride et al., 2022) or by upgrading to multi-fluid models (Meier and Shumlak, 2012;Tu and Song, 2013;Khomenko et al., 2014;Maneva et al., 2017;Popescu Braileanu and Keppens, 2022). Note that here the term ambipolar diffusion refers to the decoupling between neutral and charged particles, which causes the magnetic field to diffuse through neutral gas due to collisions between neutrals and charged particles, the latter being frozen-in into the magnetic field (Khomenko et al., 2014). In particular, Ballester et al. (2020) investigated the non-linear coupling of Alfvén and slow magneto-acoustic waves in partially ionised plasmas using a single-fluid model with ambipolar diffusion. Soler et al. (2013) and Soler, Ballester, and Zaqarashvili (2015) used two-and three-fluid models to investigate the lower solar atmosphere, showing the effects of partial ionisation and the importance of high-frequency Alfvén waves for chromospheric heating. Recently, Kuźma et al. (2020) and Pelekhata, Murawski, and Poedts (2021) studied two-fluid Alfvén waves in a partially ionised solar atmosphere. Kuźma et al. (2020) considered the case of a periodic wave driver operating at the bottom of the photosphere, and Pelekhata, Murawski, and Poedts (2021) explored impulsively generated waves in the low-atmospheric regions. They showed that these waves are able to deposit their kinetic energy in the form of heat in the chromosphere. Moreover, they showed that Alfvén waves that operated at the bottom of the photosphere do reach the chromosphere and the solar corona. In this article, we place the driver at the middle of the chromosphere and perform parametric studies for various driving amplitudes and periods. Such a position of the driver is justified as two-fluid effects are more important in the solar chromosphere, and the pure-hydrogen assumption is also more accurate in the chromosphere (Popescu Braileanu et al., 2019b). Note that a model chromosphere for studying wave heating is typically highly idealised to simplify the wave propagation and damping processes (Maneva et al., 2017;Popescu Braileanu et al., 2019b;Kuźma et al., 2020;Zhang et al., 2021;Niedziela, Murawski, and Poedts, 2021), and thus the numerical results may not quantitatively represent the processes in the solar atmosphere. The goal of this article is to examine the behaviour of Alfvén waves in a stratified two-fluid solar chromosphere that is permeated by a longitudinal magnetic field. We analyse how far upwards Alfvén waves with different periods and amplitudes can travel before they vanish and lead to a release of thermal energy and plasma outflows.
This article is organised as follows. In Section 2 the numerical model for the considered solar atmosphere is presented. The numerical results of the Alfvén waves considered in this model atmosphere are discussed in Section 3. In Section 4, we summarise our findings, formulate our conclusions, and provide a brief discussion.

Numerical Model of the Solar Atmosphere
We describe the plasma dynamics in the lower layers of the solar atmosphere by a set of twofluid equations, in which neutrals (hydrogen atoms) are described by one set of equations and ions/protons+electrons by the others (Khomenko et al., 2014;Popescu Braileanu et al., 2019a,b;Wójcik et al., 2020), neglecting the electron inertia.
Specifically, the motion of neutrals is described by the following equations: where and the motion of ions + electrons is described by where Here, γ = 5/3 is the ratio of specific heats, μ corresponds to the magnetic permeability of the medium, and the gravitational acceleration vector g = −g1 y , is oriented in the longitudinal [y] direction and has a magnitude g = 274.78 m s −2 on the Sun. Furthermore, in Equations 1-9, the subscripts i and n correspond to ions+electrons and neutrals, respectively. The symbol i,n denotes the mass densities of the respective species, while V i,n = [V i,n x , V i,n y , V i,n z ] denotes their velocities, p i and p n correspond, respectively, to the ion/proton+electron and neutral-gas pressures, B = [B x , B y , B z ] is the magnetic field, and T i,n are the temperatures, specified by the ideal-gas laws: where m i and m n denote the ion/proton+electron and neutral masses, respectively, and k B is the Boltzmann constant.
Here, the collisional momentum-[S m ] and energy-[Q i,n ] source terms are given as (Draine, 1986) and where and The coefficient of collisions between ions and neutrals is given by (e.g. Oliver et al., 2016;Ballester et al., 2018a;Soler et al., 2019, and references cited therein) where σ in corresponds to the collisional cross-section, here taken as its head-on-collision value of 1.4 × 10 −19 m 2 (Vranjes and Krstic, 2013).
The system is assumed to be invariant along the transversal [z] direction, yielding all ∂/∂z = 0, while V i,n z and B z generally differ from zero. Such a system is called quasi-1.5D and it allows Alfvén waves to propagate. Although we use a 2D mesh, we only investigate quasi-1D waves, and the horizontal span of the computational domain is small, thus reducing the computational effort.
Note that we neglected Ohmic resistance, Hall and battery effects, etc., which may lead to certain changes in the numerical results. More detailed discussions have been provided by Song and Vasyliūnas (2011) and Khomenko et al. (2014), but quantitative comparisons are not included here as they are beyond the scope of this work. We also did not include the effects of ionisation and recombination, which may affect the wave propagation and damping (Snow and Hillier, 2021;Zhang et al., 2021). Moreover, realistic ionisation and recombination are non-equilibrium processes (Leenaarts et al., 2007), and thus due to its complexity, integrating these processes is beyond the scope of our current endeavour.

Magnetohydrostatic Equilibrium
The solar chromosphere is highly dynamical (Carlsson and Stein, 1995;Wedemeyer et al., 2004). However, in this work, it is assumed that initially the solar atmosphere is in a static equilibrium, i.e. all ∂/∂t = 0 and V i = V n = 0. Note again that this assumption may lead to a further discrepancy between the present model and the realistic solar atmosphere, but it allows us to specifically study the wave-damping mechanism without being affected by other dynamical processes.
Then, from the y-components of the momentum (Equations 2 and 6), we obtain which means that gravity is balanced by the pressure gradient. Hence, with the use of the ideal-gas laws, given by Equation 10, we derive the hydrostatic gas pressures and mass densities profiles, given as Here, y r = 50 Mm is the reference level, while p 0i = 10 −2 Pa and p 0n = 3 × 10 −4 Pa are, respectively, the ion and neutral gas pressures at this level, while are the pressure scale-heights that depend on the equilibrium temperature T (y), which is specified by the semi-empirical model of Avrett and Loeser (2008).
To be clearer, several important quantities are shown in Figure 1. First, the temperature is shown in Figure 1 (top). Note that T (y = 0) ≈ 5800 K. Higher up the temperature declines with height [y] above the solar surface, reaching its minimum of 4341 K at y ≈ 0.6 Mm. Then, T grows with y to about 6 × 10 5 K at y = 5 Mm. Equilibrium ion (solid line) and neutral (dashed line) mass densities are displayed in Figure 1 (upper-middle). Note that there are about 10 3 times more neutrals than ions at the bottom of the photosphere, y = 0 Mm. Higher up, the relative number of neutrals in comparison to ions declines with height y, and the solar corona consists essentially of ions. The hydrostatic equilibrium is overlaid by a uniform (i.e. force-and current-free) longitudinal magnetic field: B = B y1y . We consider the case of B = 10 G. The corresponding total Alfvén speed c A = B y / μ( i + n ) is displayed in Figure 1 (lower-middle). Note that c A = 0.1 km s −1 at y = 0 Mm and c A = 3.5 km s −1 at y = 1 Mm. Higher up, it grows with height y, attaining its value of 540 km s −1 at y = 5 Mm. Figure 1 (bottom) displays neutral-ion collision time [τ ni = n /α in ] vs height. Note that τ ni (y = 0) ≈ 2 × 10 −6 s. It grows dramatically with y such that at y = 2 Mm τ in is about 2 × 10 −3 s. Above the transition region, τ in increases gently with height until at y = 30 Mm τ in ≈ 0.1 s. Similar values of the neutral-ion collision time were obtained by Alharbi et al. (2022) based on the VAL III atmospheric model.
In this work, neglecting ionisation and recombination results in a more ideal background field for investigating the wave-propagation process, although longitudinal flows may change the ionisation fraction distribution to some extent. Nevertheless, with a simplified model similar to that of Zhang et al. (2021), the effects of including ionisation and recombination are being studied.

Monochromatic Driver
The above-described equilibrium is perturbed by a mono-periodic driver specified in the boundary values of the z-components of ion and neutral velocities at height y = y 0 . Although a broadband perturbation may reproduce a more realistic energy source (Tu and Song, 2013), the present driver allows us to specifically investigate the energy-transport process with several given wave periods. This driver is set in the form of a plane-wave front, given by where V 0 is the amplitude of the driver and P d stands for its period. Below, we set and hold fixed y 0 = 1 Mm and consider several values of P d , from 2.5 s to 10 s, and V 0 taking either 0.1 km s −1 or 1 km s −1 . The amplitudes adopted here approximately resemble the transverse motions observed in the photosphere (Matsumoto and Shibata, 2010) and the chromosphere (De Pontieu et al., 2007). The driver does not directly change the other velocity components, nor the other quantities. While the resulting perturbations have non-vanishing components V i z and B z , they correspond to Alfvén waves. Moreover, V n z is the transversal component of the neutral velocity. This transverse motion of the neutrals is caused by the collisional interactions between neutrals and charged particles, as only the charged particles interact with the magnetic field. Note again that in this work, we use an idealised model solar atmosphere and a simplified two-fluid model, and more detailed investigations and improvements are relegated to the future.

Dispersion Relation
The attenuation of the Alfvén waves is our main interest in this article. To measure it we define the attenuation length [λ d ] as For comparison, we can first calculate λ d from the dispersion relation for linear Alfvén waves propagating in an ion-neutral two-fluid homogeneous medium (Zaqarashvili, Khodachenko, and Rucker, 2011;Ballester et al., 2018a) where Here, i 2 = −1 and k = k r + ik i is the complex wave number from which we can calculate the wavelength [λ = 2π/k r ] and the attenuation length [λ d = 2π/k i ] of the excited Alfvén waves. The dispersion relation used in this section is not strictly applicable to the atmosphere, which is gravitationally stratified, but under local conditions it can be used as a valid first-order approximation for linear waves. Also, note that there are more sophisticated descriptions for an electron-ion-neutral three-fluid model (Song, Vasyliūnas, and Ma, 2005) or an ion-neutral multi-species model (Zaqarashvili, Khodachenko, and Rucker, 2011), but they are not further discussed for simplicity. Figure 2 shows the analytically obtained values of the wavelength [λ] (top) and the attenuation length [λ d ], (bottom). For instance, at y = 1 Mm, for the mass densities we took i = 7.9 × 10 −9 kg m −3 , n = 4.7 × 10 −9 kg m −3 . Then, we infer that at this height for P d = 5 s, λ ≈ 14 km and λ d = 5.3 × 10 5 km. The analytical predictions for linear Alfvén waves suggest that these waves are weakly damped by ion-neutral collisions. As a result, we expect that the linear Alfvén waves will only be weakly damped in the chromosphere. Numerical results for non-linear Alfvén waves are also discussed in detail in the next section.

Numerical Results
In this section, we present the results of 1.5-D numerical simulations performed with the use of the JOANNA code (Wójcik et al., 2019, which numerically solves the initialboundary value problem for the two-fluid equations. A second-order accurate spatial reconstruction (Toro, Hidalgo, and Dumbser, 2009) and a third-order accurate Super Stability Preserving Runge-Kutta method (Durran, 2010) are used. Numerical fluxes are given by using the HLLD approximate Riemann solver (Miyoshi and Kusano, 2005), and the divergence of the magnetic field is controlled by the cleaning method of Dedner et al. (2002). The CFL number (Courant, Friedrichs, and Lewy, 1928) is set to 0.9 to determine the time steps.
While the initial field in a magneto-hydrostatic equilibrium at t = 0 s has been introduced above, we now further introduce the discretisation. For the wave period P d = 2.5 s, P d = 5 s, and P d = 10 s we set a uniform grid with a cell size of y = 0.5 km in the longitudinal direction, which extends over y 0 ≤ y ≤ 2 Mm. We performed some simulations for the  spatial resolution with the grid size y = 0.1 km, and for the current wave frequencies the numerical results do not differ much from those obtained for y = 0.5 km. For larger values of y the computational box was covered by 256 grid points up to y = 60 Mm. At the bottom and top boundaries, we set and hold fixed in time all plasma quantities to their equilibrium values. Additionally, the bottom boundary is overlaid by the driver given by Equation 20.

Numerical Disturbances
A standard polynomial reconstruction scheme tends to produce spurious disturbances in gravitationally stratified fluid flows, where the density and pressure are exponentially decreasing with altitude, and the disturbances may affect the numerical results (Käppeli and Mishra, 2016;Krause, 2019;Brchnelova et al., 2022;Käppeli, 2022). Therefore, we first examine the numerically produced disturbances.
Guided by this idea, we ran a simulation without any driver and for the present grid resolution y = 0.5 km. Figure 3 presents the longitudinal component of the ion velocity (top) and the relative perturbed ion temperature (bottom), where T = T i (y, t = 0 s) = T n (y, t = 0 s) is taken from the temperature of the average quiet-Sun chromosphere model of Avrett and Loeser (2008). Note that max(V iy ) ≈ 0.015 km s −1 occurs in the solar corona. This value is much less than the local Alfvén speed, which reaches about 540 km s −1 at y = 5 Mm (Figure 1, lower-middle). Similarly, the max( T i /T ) ≈ 0.01% is sufficiently low (Figure 3, bottom). More specifically, in both the longitudinal velocity and temperature, we can see some disturbances originating from the transition region, where the plasma quantities change by several orders of magnitude. However, the disturbances are not dominant and continue to decrease over time. As a result, we conclude that for the chosen grid resolution the numerical disturbances are negligibly small, and thus the following simulations can start at t = 0 s.

Small-Amplitude Alfvén Waves
In this section we consider small-amplitude Alfvén waves. Figure 4 displays time-distance plots for V iz in the case of V 0 = 0.1 km s −1 , which is linear as its amplitude is significantly smaller in comparison with the local Alfvén speed. Both for P d = 2.5 s (top) and P d = 5 s (middle) Alfvén waves are more strongly damped than for P d = 10 s (bottom) but they all reach the solar corona. In general, the wave damping in these linear cases is small, and thus the numerical results agree with the analytical prediction from the dispersion relation of Equation 22, which suggests that linear waves are only weakly damped. This agreement further verifies the code being used. Moreover, note that for all considered wave periods, V i contours show clear patterns of standing waves, which result from the reflection of Alfvén waves at the transition region. The reflection of Alfvén waves has been well explained by Tu and Song (2013). Although only part of the energy is reflected back into the chromosphere, it still enhances the waves, which will be more evident in the non-linear cases.

Large-Amplitude Alfvén Waves
In the previous subsection, we have seen that the small-amplitude Alfvén waves are weakly damped. Martínez-Gómez, Soler, and Terradas (2018) discussed non-linear Alfvén waves propagating in homogeneous partially ionised media, and showed that damping may also relate to the wave amplitude. Therefore, in the present section, we consider large-amplitude waves. Specifically, we consider the case of the driver amplitude V 0 = 1 km s −1 . In the following simulations, the transverse velocity V i z is not a small value in comparison with the Alfvén speed. For instance, the maximum value of ratio V i z /c A is larger than 20% for P d = 10 s (not shown). Therefore, the linear theory may no longer be accurate and extra damping may be produced due to the non-linearity. Non-linear Alfvén waves can generate steep wave fronts and drive longitudinal compressive perturbations (Montgomery, 1959;Hollweg, 1971;Nakariakov, Ofman, and Arber, 2000;Vasheghani Farahani et al., 2012). As a result, the non-linear Alfvén waves may experience direct damping at steep wave fronts and indirect damping by the driven magneto-acoustic waves, which are affected by collisional damping and shock damping (Popescu Braileanu et al., 2019b;Zhang et al., 2021).
We can see the non-linear behaviours of the P d = 10 s case clearly in Figure 5. At the location of the driver (y = 1 Mm) the amplitude of the magneto-acoustic waves is about 0.1 km s −1 . Higher up, approximately when the Alfvén speed surpasses the sound speed, the longitudinal velocity starts to significantly increase; at about y = 1.4 Mm the amplitude has grown about fivefold; in the meantime, shocks also appear and the amplitude stops increasing. Since the mass density is gravitationally stratified, the energy carried by the magnetoacoustic waves in fact decreases when the mass density decreases and the velocity does not increase. Therefore, a part of the wave energy has to be deposited in the plasma. Moreover, the strong longitudinal wave steepens to shocks, which may cause sudden changes of the density, which then affect the Alfvén speed, and this partly explains that the interaction between different (Alfvén and acoustic) wave modes causes relatively irregular wave patterns. In contrast, when the period is short (P d = 2.5 s), the frequency of the magneto-acoustic wave is close to 1 Hz, and it experiences stronger non-linear damping even before shocks appear (Popescu Braileanu et al., 2019b). Therefore, in the P d = 2.5 s case, both the Alfvén and acoustic waves are regular and smooth. Note that in the transition region, the amplitudes of both the Alfvén waves and magneto-acoustic waves barely amplify, because part of their energy is reflected, which will be further discussed below.
In the following paragraphs, we do not discuss much about the details of the damping mechanism of the non-linear waves, as non-linear (M)HD waves and shocks were already discussed extensively (Ballester et al., 2020;Snow and Hillier, 2020;Zhang et al., 2021;Snow and Hillier, 2021). Instead, we focus on the impact resulting from these waves.

Longitudinal Flows
It is known that the longitudinal magnetic pressure gradient [−∂(B 2 z /2μ)/∂y] enables Alfvén waves to drive a longitudinal flow V i y (Hollweg, 1971). In a low plasma-β region V iy corresponds to (slow) magneto-acoustic waves, which propagate approximately at c s , and V i z is always associated with Alfvén waves, moving at c A . When the Alfvén waves are generated by a small-amplitude driver, e.g. with V 0 = 0.1 km s −1 as in the above simulations, the longitudinal flow is certainly small (not shown). However, for larger driver amplitudes, e.g. for V 0 = 1 km s −1 as considered here, the longitudinal flows ( Figure 6) are not negligible. In particular, in the case of a wave period P d = 10 s the longitudinal flow attains its maximum value of about 3.3 km s −1 at the top of the chromosphere (Figure 6, bottom). Note that in Figure 6, there are some initial disturbances rising from the bottom boundary, which are the longitudinal flows initially driven from the bottom of the domain, before the waves become stationary. These disturbances become weaker over time and eventually, the wave patterns become regular and stationary. Moreover, we can also see reflections in the longitudinal-velocity contours. The reflection is stronger in the higher-frequency case (P d = 2.5 s), and thus the amplitude of the longitudinal velocity largely decreases in the transition region (which is counter-intuitive since there is a strong density stratification). Moreover, Figure 7 illustrates the longitudinal component of the ion velocity averaged over time [t ]: where t p is a duration of the transient phase, chosen as 200 s and t max = 500 s. For all the presented cases of P d and in the range 1 Mm y 0 < 2 Mm, V i y t is positive, which corresponds to plasma outflows. The maximum values reach 0.1 km s −1 at y ≈ 2 Mm for P d = 10 s. Such plasma outflows follow the general trend noted also in recent observational data (Kayshap, Banerjee, and Srivastava, 2015). Note that for P d = 2.5 s, V i y t attains its lowest values of about 0.04 km s −1 , which means that a larger value of P d results in a larger net outflow. It is worth noting that for P d = 2.5 s, the longitudinal flow velocity is higher than in the other two cases in the chromosphere. This is probably because this shorter wave experiences a stronger reflection in the transition region and stronger non-linear damping in the chromosphere (Popescu Braileanu et al., 2019b). With a longer wave period, the waves have larger amplitudes above the chromosphere, which is why they lead to stronger outflows above the chromosphere.

Plasma Heating
We now investigate plasma heating. In Figure 8 (top) we can see that the wave pattern is similar to that of the small-amplitude waves. However, it is more irregular here due to the non-linear behaviours mentioned above. Not surprisingly, the absolute value of the longitudinal velocity drift between ions and neutrals [|V i y − V n y | (Figure 8, upper-middle)] attains its largest values in the higher layers of the chromosphere, as the coupling becomes weaker with (much) lower density, but it is about one order of magnitude smaller than the absolute value of the transverse velocity drift (Figure 8, lower-middle). As a result, ion-neutral collisions produce less thermal energy from slow magneto-acoustic waves than from Alfvén waves. However, as the acoustic waves steepen to shocks, they contribute to the overall heating by non-linear damping (Ballester et al., 2018b;Popescu Braileanu et al., 2019b). Note that in the longitudinal-velocity drift contour we can see some disturbances coming down from the transition region. The disturbances are much smaller than the waves, as already explained above in Section 4.1. T i /T reaches its maximum value of about 12% near y = 1.35 Mm (Figure 8, bottom), which is approximately the location where the amplitude of the magneto-acoustic wave stops increasing and even decreases ( Figure 5, top). This finding also suggests that the non-linear damping of the magneto-acoustic wave is the dominant heating mechanism, at least in this case. The collisional heating term [Q v ] is displayed in Figure 9. With the large amplitude, the collisional heating is stronger, although it is still too small to be sufficient for the observed Figure 8 Time-distance plots for V i z top, |V i y − V n y | upper-middle, |V i z − V n z | lower-middle, and relative perturbed ion temperature T i /T bottom for V 0 = 1 km s −1 and P d = 10 s. temperature increases. The reason for this is that the periods considered here are too large in comparison with the collision time. This was well explained by Song, Vasyliūnas, and Ma (2005), who showed that the collisional damping is the strongest in the collision-time range between the ion-neutral collision time and neutral-ion collision time.
To further quantify the heating, we calculate the time average of the relative perturbed temperature [ T i /T t ] given by for all heights y from y 0 = 1 Mm to the upper chromosphere at y = 1.9 Mm, and for the three considered values of P d . The results are shown in Figure 10, which illustrates T i /T t in the cases of P d = 2.5 s, P d = 5 s, and P d = 10 s, for V 0 = 1 km s −1 . It is worth noting that stronger heating takes place for waves with P d = 2.5 s. This may be related to the stronger (acoustic wave) reflection that keeps more energy in the chromosphere and the stronger non-linear damping of the higher-frequency waves. Figure 11 illustrates the relative perturbed temperature [ T i /T t,y ] averaged over all simulation time and height from y = y 0 to y = y 0 + λ H obtained from where λ H is the heating length, which we fix as λ H = 2 Mm. Note that heating of the solar atmosphere by the damped Alfvén waves falls off with P d , which agrees with the analytical calculations presented by Song, Vasyliūnas, and Ma (2005), Vranjes et al. (2008), Zaqarashvili, Khodachenko, and Soler (2013). These results in fact also support the conclusion Figure 11 The relative temperature averaged over time [ T i /T t,y ] from t = 200 s to t = 500 s and height from y = y 0 to y = y 0 + λ H , versus P d for V 0 = 1 km s −1 .
of De Pontieu et al. (2007), who suggested that the strong Alfvén waves in the chromosphere could provide sufficient energy to heat the quiet solar corona and to accelerate the solar wind.

Summary and Conclusions
We constructed an idealised magnetohydrostatic model of the solar atmosphere in which a monochromatic driver, while acting in the middle chromosphere, excites Alfvén waves. This atmosphere is described by a set of simplified two-fluid equations for ions + electrons and neutrals, treated as two separate fluids that interact with each other by collisions, while omitting the Hall and battery effects, ionisation/recombination, etc. These two-fluid equations are solved numerically with the use of the JOANNA code (Wójcik et al., 2019).
We found that small-amplitude Alfvén waves are very weakly damped in the chromosphere, whereas, with a stronger velocity-amplitude driver, the driven large-amplitude Alfvén waves behave differently. Mainly, the large-amplitude non-linear waves experience stronger damping than the small-amplitude waves, resulting in significant energy thermalisation in the chromosphere. Specifically, non-linear Alfvén waves exhibit steep wave fronts and they drive magneto-acoustic waves. As a result, these Alfvén waves experience direct non-linear damping, and the magnetic-acoustic waves also experience shock damping. Moreover, the large-amplitude Alfvén waves drive stronger longitudinal flows toward the corona. Such flows were recently detected by Kayshap, Banerjee, and Srivastava (2015) who reported plasma outflows of about 2 km s −1 above the transition region. Our results indicate that the large-amplitude Alfvén waves may be the mechanism driving the solar wind. Data Availability All data generated or analysed during this study are included in this published article.

Competing interests The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.