Multifluid approach to highfrequency waves in plasmas. III. Nonlinear regime and plasma heating
Abstract
The multifluid modelling of highfrequency waves in partially ionized plasmas has shown that the behavior of magnetohydrodynamics waves in the linear regime is heavily influenced by the collisional interaction between the different species that form the plasma. Here, we go beyond linear theory and study largeamplitude waves in partially ionized plasmas using a nonlinear multifluid code. It is known that in fully ionized plasmas, nonlinear Alfvén waves generate density and pressure perturbations. Those nonlinear effects are more pronounced for standing oscillations than for propagating waves. By means of numerical simulations and analytical approximations, we examine how the collisional interaction between ions and neutrals affects the nonlinear evolution. The friction due to collisions dissipates a fraction of the wave energy, which is transformed into heat and consequently rises the temperature of the plasma. As an application, we investigate frictional heating in a plasma with physical conditions akin to those in a solar quiescent prominence.
1 Introduction
Alfvén waves are usually classified in two categories depending on their velocity amplitudes. Smallamplitude waves have velocity amplitudes that are much smaller than the Alfvén speed. Conversely, the velocity amplitudes of largeamplitude waves are not negligible in comparison with the Alfvén speed. In MartínezGómez et al. (2016) and MartínezGómez et al. (2017) (hereafter, Papers I and II, respectively), we studied the properties of smallamplitude waves in fully and partially ionized plasmas of the solar atmosphere by considering the linear regime of a multifluid model. The goal of those works was to investigate the effects of the collisional interactions between the different species in multicomponent plasmas. In the present paper, we extend those studies by incorporating the nonlinear effects that arise when largeamplitude perturbations are considered.
Solar prominences and the solar wind are examples of solar plasmas in which largeamplitude waves have been detected. Largeamplitude oscillations of solar prominences are typically triggered by nearby flares (Moreton & Ramsey, 1960; Ramsey & Smith, 1966): these events produce Moreton and EIT waves (Moses et al., 1997; Thompson et al., 1998) that impact on the prominence causing the whole structure to vibrate during a few periods. These largeamplitude oscillations are rare events but in the last years a growing number of observations has been reported (see, e.g., Eto et al., 2002; Jing et al., 2003; Okamoto et al., 2004; Gilbert et al., 2008; Luna & Karpen, 2012). In addition, Coleman (1968), Belcher et al. (1969) and Belcher & Davis (1971) detected largeamplitude waves in the solar wind. Since those observations, the properties of nonlinear magnetohydrodynamic (MHD) waves have been extensively studied by, e.g., Hollweg (1971), Cohen & Kulsrud (1974), Lou (1993) or Ruderman (2006), and their role in processes like the acceleration of the solar wind or the heating of the solar atmosphere has been investigated by, e.g., Barnes & Hollweg (1974), Esser et al. (1986), Ofman & Davila (1998), Voitenko & Goossens (2002) or Suzuki (2008), among many others.
The theoretical study of nonlinear waves is more involved than that of their linear counterparts and is typically performed by means of numerical simulations (see, e,g., Murawski & Roberts, 1993; Oliver et al., 1998). Some more recent numerical results can be found in Matsumoto & Shibata (2010), who studied Alfvén waves driven by photospheric motions, Suzuki (2011), who investigated solar and stellar winds driven by Alfvén waves, or Karpen et al. (2017), whose results suggest that coronalhole jets are a possible origin for nonlinear Alfvén waves in the interplanetary medium.
Nevertheless, analytical results can also be obtained if certain approximations are taken. A common analytical procedure is to assume a perturbative expansion, where the variables that describe the properties of the plasma are expressed as a sum of a background value plus a series of terms that represent the linear and higherorder perturbations. The series is truncated at some given order and systems of equations are derived for the perturbations, while higherorder effects are neglected. This procedure was followed, e.g., by Hollweg (1971), who studied secondorder effects of Alfvén waves, or by Rankin et al. (1994), Tikhonchuk et al. (1995) and Verwichte et al. (1999), who examined the properties of up to thirdorder perturbations. Those works have shown that nonlinear Alfvén waves induce a ponderomotive force that causes variations in the density and pressure of the plasma, in contrast with the incompressibility of linear Alfvén waves. In addition, thirdorder effects also produce a steepening of the wave and the generation of higher harmonics.
In most of the works mentioned in the previous paragraphs the plasma is considered to be fully ionized and treated as a singlefluid. The assumption of fully ionization is valid for the solar corona and the solar wind, where the presence of neutral particles is negligible. However, it is not applicable to other regions of the solar atmosphere, such as the chromosphere or prominences, in which neutrals are the dominant component of the plasma and have a dramatic effect on the properties of MHD waves (see, e.g., Piddington, 1956; Watanabe, 1961; Haerendel, 1992; Soler et al., 2013a). In addition, as shown in Papers I and II, the use of singlefluid models is only appropriate when the phenomena of interest is associated with low frequencies, i.e., much lower than the ion cyclotron frequencies in fully ionized plasmas or the ionneutral collisions frequencies in partially ionized plasmas. Conversely, at higher frequencies, multifluid approaches are required due to the fact that the components of the plasma are not strongly coupled and react to perturbations in different timescales.
In the present work, the multifluid model described in Paper I is applied to the investigation of nonlinear waves in partially ionized plasmas, paying special attention to the heating due to ionneutral collisions. The issue of heating is of great interest in the solar atmosphere (see, e.g., Goodman, 2011; Parnell & De Moortel, 2012; Song & Vasyliūnas, 2011; Khomenko & Collados, 2012; Tu & Song, 2013; Gilbert, 2015; Heinzel, 2015; Arber et al., 2016; Soler et al., 2016, 2017). It has been shown that Alfvénic waves can transport a huge amount of energy from the photosphere to higher layers of the solar atmosphere (Tomczyk et al., 2007; McIntosh et al., 2011; Srivastava et al., 2017). However, it remains unclear whether all the energy carried by the waves is deposited in the plasma. A dissipative mechanism is required to transform that energy into heat and, in the case of partially ionized plasma, the ionneutral collisional interaction is one of the possible mechanisms. The topic of heating by means of ionneutral collisions was briefly examined in Paper II when smallamplitude perturbations were studied. However, since heating is a nonlinear effect with a quadratic dependence on the velocity drifts, as shown by Equation (1) of Paper II, it is expected to have a more relevant role when largeamplitude waves are considered. The model used here also considers Coulomb collisions, magnetic diffusivity and the effects of Hall’s current, which accounts for the cyclotron motion of ions and has been shown to be of great relevance in weakly ionized plasmas (see, e.g., Pandey & Wardle, 2008).
The outline of this paper is as follows. In Section 2, the effect of partial ionization on nonlinear standing waves is investigated: numerical simulations are performed for the case of a plasma with prominence conditions. In Section 3, largeamplitude impulsive perturbations are considered and the heating due to ionneutral collisions is examined. In Section 4, we study the propagation of nonlinear Alfvén waves generated by a periodic driver. Finally, Section 5 summarizes the results of this work. As complementary content, the appendix includes analytical results for the case of partially ionized twofluid plasmas.
2 Standing waves
In this section, nonlinear standing waves in a uniform and static partially ionized plasma are analyzed. The temporal evolution is governed by the equations detailed in Section 2 of Paper I. In short, we use a combination of the 5moment transport equations (Schunk, 1977) for each species of the plasma, the induction equation obtained from Faraday’s law, and a generalized Ohm’s law that includes Hall’s term, the Biermann battery term (related to the gradient of electronic pressure) and magnetic resistivity or Ohm’s diffusion. Interested readers are referred to Paper I for the detailed derivation and discussion of the governing multifluid equations. Due to the complexity of the equations, we perform 1.5D numerical simulations with a modified version of the MolMHD code (Bona et al., 2009).
Figure 1 shows the results of a simulation in a plasma with conditions that correspond to a quiescent prominence core at an altitude of 10,000\ \mathrm{km} over the photosphere and with gas pressure of P_{g}=0.005\ \mathrm{Pa}, according to Heinzel et al. (2015). The plasma in the cool prominence is composed of three different species: protons, neutral hydrogen and neutral helium, denoted by the subscripts p, \mathrm{H}, and \mathrm{He}, respectively. Ionized helium has a residual abundance at the cool temperatures of prominence cores. So, for simplicity, we assume helium to be fully neutral. The presence of ionized helium is, however, important for prominencetocorona transition region temperatures. At t=0 the number densities of protons, neutral hydrogen, and neutral helium, and the temperature are uniform. Their values are n_{p}=1.4\times 10^{16}\ \mathrm{m^{3}}, n_{\mathrm{H}}=2\times 10^{16}\ \mathrm{m^{3}}, n_{\mathrm{He}}=2\times 10^{15}\ \mathrm{m^{3}}, and T_{0}=10,000\ \mathrm{K}, respectively. The collision frequencies of protons with neutral hydrogen, of protons with neutral helium, and of hydrogen with helium are \nu_{p\mathrm{H}}\approx 270\ \mathrm{Hz}, \nu_{p\mathrm{He}}\approx 3.5\ \mathrm{Hz}, and \nu_{\mathrm{HHe}}\approx 5.2\ \mathrm{Hz}, respectively. These values are computed from the friction coefficients given by Equation (4) of Paper II. The collision frequency, \nu_{st}, between two species s and t and the friction coefficient, \alpha_{st}, are related by \nu_{st}=\alpha_{st}/\rho_{s}, where \rho_{s}=m_{s}n_{s} is the density of the species s and m_{s} is the particle mass. A uniform background magnetic field, \bm{B_{0}}, along the xdirection is considered. A typical value of the magnetic field strength in quiescent prominences is B_{0}=10\ \mathrm{G}. The fundamental standing mode of the transverse Alfvén waves is excited by applying the initial perturbation
V_{s,y}(x,t=0)=V_{y,0}\cos(k_{x}x),  (1) 
to every species s of the plasma, where V_{s,y} is the ycomponent of the velocity and k_{x} is the longitudinal wavenumber. No initial perturbation is applied to the other variables, i.e., V_{s,x}(x,t=0)=V_{s,z}(x,t=0)=0,\bm{B}(x,t=0)=\bm{0}, and \rho_{s}(x,t=0)=P_{s}(x,t=0)=0. These initial conditions generate circularly polarized Alfvén waves. The domain we have chosen for this simulation is x\in[l,l], with l=2.5\times 10^{5}\ \mathrm{m}. Thus, the wavenumber of the fundamental mode is k_{x}=\pi/(2l)=\pi/(5\times 10^{5})\ \mathrm{m^{1}} and the boundary conditions impose that the three components of the velocity are equal to zero at x=\pm l, while the rest of variables are extrapolated. The amplitude of the perturbation is given by V_{y,0}=2.5\times 10^{2}c_{\mathrm{A}}, where c_{\mathrm{A}}=\bm{B_{0}}/\sqrt{\mu_{0}\rho_{p}} is the Alfvén speed, with \mu_{0} the vacuum magnetic permeability. Note that the present definition of Alfvén speed only takes into account the density of ions. For the parameters given above, the Alfvén speed is c_{\mathrm{A}}\approx 184\ \mathrm{km\ s^{1}}.
The top row of Figure 1 shows several snapshots of the evolution of the ycomponent of the velocity, which is perpendicular to the background magnetic field. Although they have not been represented here, the initial condition given by Equation (1) also generates perturbations on the ycomponent of the magnetic field and the zcomponents of both the velocity and the magnetic field. This means that, as time advances, the oscillation plane of the waves rotates. However, for the range of frequencies studied here, the rotation is very slow and the amplitudes of V_{s,z} and B_{z} are much smaller than those of V_{s,y} and B_{y} throughout the whole duration of the simulations. In addition, B_{y} has a temporal and spatial phase shift with respect to V_{y}, as expected for a standing Alfvén wave, but it has the same behavior in terms of frequency and damping of its oscillations. For those reasons, to analyze the properties of the transverse waves we focus only on the ycomponent of the velocity.
Initially, the three species of the plasma have the same velocity but, since the coupling between them is not perfect, some small velocity drifts appear when the Alfvén wave starts its oscillation. As time advances, the collisional friction causes the damping of the wave, as is better illustrated by the top left panel of Figure 2. This is the same behavior as the one already explained in Paper II for smallamplitude waves. Nevertheless, due to the much larger amplitude of the perturbation used in the present investigation, the nonlinearities are not negligible and perturbations along the direction of the background magnetic field are also excited. Thus, the second row of Figure 1 displays the xcomponent of the velocity, normalized with respect to the amplitude of the driver, V_{y,0}, at various time steps. Apart from the smaller amplitude of V_{x} in comparison with V_{y}, the main difference is that its wavenumber is twice the wavenumber of the initial perturbation and there is a spatial phase shift: while V_{y} is proportional to \cos(k_{x}x), V_{x} is proportional to \sin(2k_{x}x). Furthermore, the right top panel of Figure 2 shows that the oscillation in V_{x} does not attenuate as fast as the oscillation in V_{y}.
The third and bottom rows of Figure 1 show the relative variation of density, defined as the ratio between the perturbation in density and its background value, i.e., \Delta\rho/\rho_{0}, with \Delta\rho\equiv\rho(x,t)\rho_{0}, and the ratio between the perturbation of temperature and the initial temperature, \Delta T/T_{0}, with \Delta T\equiv T(x,t)T_{0}, respectively. The temperature of each species s is computed from its pressure and number density through the ideal gas law, P_{s}=n_{s}k_{\mathrm{B}}T_{s}, where k_{\mathrm{B}} is Boltzmann’s constant. It can be checked that both \Delta\rho/\rho_{0} and \Delta T/T_{0} are proportional to \cos(2k_{x}x). The relative variation of density shows that matter accumulates at the center of the domain and is displaced from the ends during the first steps of the simulation but later this process is reversed and a oscillation appears. The amplitude of this variation of density is around 2\% of the initial background value. Finally, the bottom panels show that the average temperature of the plasma (denoted by a horizontal dotted line) rises as time advances. The main reason for this increase is the dissipation of the kinetic energy of the initial perturbation which is transformed into heat by means of ionneutral collisions.
More details of the simulation can be analyzed by inspecting Figure 2, where the temporal evolution of the same variables displayed in Figure 1 at selected representative points of the domain is plotted. The representative point for V_{x} is different from the position chosen for the rest of variables because x=0 is a node for this component of the velocity. Hence, a better location to analyze V_{x} is x=l/2.
By fitting the oscillation displayed in Figure 2(a) with an exponentially decaying sinusoidal of frequency \omega, we find that \omega\approx 0.67\ \mathrm{rad\ s^{1}} (which corresponds to a period of 9.4\ \mathrm{s}). This frequency agrees with the result obtained by solving the dispersion relation derived in Paper II for linear perturbations, namely Equation (16) of Paper II. If we compare the collision frequencies between the different fluids with the oscillation frequency divided by 2\pi (it is common practice to compare directly \omega with \nu_{st} but, rigorously speaking, this comparison is not correct because those quantities are expressed in different units), we find that \omega/(2\pi)<\nu_{st}. This fact explains why the three species oscillate with almost the same velocity but there is still some damping due to friction.
The top right panel of Figure 2 shows a wave in the xcomponent of the velocity that seems to be composed of at least two different oscillation modes. The longitudinal motion is dominated by a mode that oscillates with a frequency much lower than the frequency of the transverse oscillation and is weakly damped. The analysis of the longitudinal oscillation reveals that the frequencies of the two modes are \omega_{1}\approx 0.16\ \mathrm{rad\ s^{1}} and \omega_{2}\approx 1.34\ \mathrm{rad\ s^{1}}. We show later that these frequencies are related to the weighted mean sound speed of the whole fluid and the Alfvén speed (modified by the inclusion of the density of neutrals), respectively.
Panel (c) of Figure 2 shows the temporal evolution of \Delta\rho/\rho_{0} at x=0. It can be seen that at the central point of the simulation domain the density rises during the initial seconds, it reaches a maximum and then the fluctuation can be described as the composition of an oscillation and a decreasing trend with time. It can be checked that the frequency of the density oscillation coincides with that of the dominant mode in V_{x} and that there is a slight temporal phase shift between \Delta\rho/\rho_{0} and V_{x}. Finally, panel (d) shows a growing trend of the temperature at x=0, combined with an oscillation similar to that found in the density. This increase of temperature is a consequence of the friction due to ionneutral collisions. A fraction of the energy of the Alfvén wave is transformed into heat and, thus, the internal energy of the plasma grows.
The previous results have been obtained for a case with a strong coupling between the three fluids of the plasma. It is interesting to repeat the simulations but when the interaction between fluids is weaker. This can be achieved by considering a wave with \omega/(2\pi)>\nu_{p\mathrm{He}}. To that goal, we perform a simulation with a larger wavenumber, k_{x}=\pi/(5\times 10^{3})\ \mathrm{m^{1}}. The dispersion relation (see Paper II) predicts a frequency \omega\approx 75.14\ \mathrm{rad\ s^{1}} for the Alfvén wave, which is higher than 2\pi\nu_{p\mathrm{He}} and 2\pi\nu_{\mathrm{HHe}}, but lower than 2\pi\nu_{p\mathrm{H}}. The results of this simulation are displayed in Figure 3. Remarkable differences with respect to the previous case are found. Now, the attenuation is stronger than in Figure 2 and the Alfvén wave dissipates quickly. Moreover, neutral helium is found to be decoupled from the other species. In contrast with the previous case, panel (b) shows that the wave in the xcomponent of velocity is more attenuated with time and, in addition, only one oscillation mode can be clearly noticed instead of the two modes present in the first simulation. Some hints of the second mode may be found during the first instants of the motion but it disappears fast. Again, it is evident that neutral helium is not as strongly coupled to protons and neutral hydrogen as before. The decoupling of neutral helium from the rest of species is a purely multifluid effect that cannot be captured with the usual singlefluid treatments.
Figure 3(c) shows that the density only increases at the center of the domain during a very short time. Then, the relative variation of density becomes negative and oscillates about \Delta\rho/\rho_{0}\approx0.02. Hence, the net result of this nonlinear effect is that matter is displaced from the central part of the domain and directed towards the ends. This behavior may be related to the increase of the fluid pressure which is associated with the initial fast grow of temperature shown in panel (d). The quick rise of temperature and pressure is caused by the fast dissipation of the Alfvén wave due the collisional friction. This issue is addressed with more detail later. It is also remarkable that during the first steps of the simulation, neutral helium reaches a larger temperature than the other two fluids and then tends to a thermal equilibrium with them. This is all caused by collisions, which tend to equalize the temperatures of all species in a timescale of the order of the collision frequency (Spitzer, 1956).
To gain a better understanding of the nonlinear effects presented up to this point, it would be useful to derive some analytical expressions from the multifluid equations. However, a threefluid system is quite complex for this goal and it would be difficult to extract some clear conclusions. A simpler scenario that can be investigated analytically is the case of partially ionized plasmas composed of only two distinct fluids. A detailed analysis of such simpler scenario can be found in Appendix A. Here, we just mention its main conclusions:

a standing Alfvén wave nonlinearly generates two secondorder perturbations in density, as well as in pressure and the longitudinal component of the velocity. The wavenumber of those perturbations, \kappa, is twice the wavenumber of the original wave, i.e., \kappa=2k_{x};

the frequencies of the secondorder perturbations are given by 2\widetilde{c}_{S}k_{x} and 2\widetilde{c}_{\mathrm{A}}k_{x}, where \widetilde{c}_{S} and \widetilde{c}_{\mathrm{A}} are the effective sound speed and the modified Alfvén speed, respectively. The former is given by the positive square root of
\widetilde{c}_{S}^{2}=\frac{\sum_{t}\rho_{t}c_{S,t}^{2}}{\sum_{t}\rho_{t}}, (2) where c_{S,t}=\sqrt{\gamma P_{t,0}/\rho_{t,0}} is the sound speed of species t, with \gamma the adiabatic constant and P_{t,0} the equilibrium value of pressure. The modified Alfvén speed is given by
\widetilde{c}_{A}=\frac{c_{\mathrm{A}}}{\sqrt{1+\sum_{t\neq p}\chi_{t}}}, (3) where \chi_{t}=\rho_{t}/\rho_{p};

and, if \widetilde{c}_{S}^{2}\ll\widetilde{c}_{\mathrm{A}}^{2}, the relative variation of density is dominated by the mode associated with the effective sound speed and is proportional to V_{y,0}^{2}/\widetilde{c}_{S}^{2}.
These conclusions are in good agreement with the numerical results represented in Figures 13. However, it must be noted that they correspond to a case in which the coupling between all the species of the plasma is strong. We are also interested in the study of scenarios with weaker couplings. Hence, we next perform a series of simulations to investigate the dependence of the secondorder perturbations on the wavenumber. The results of this study are represented in Figure 4, where the normalized oscillation frequency of waves, \omega_{R}/\omega_{\mathrm{A}}, is plotted as a function of the wavenumber k_{x}. The normalization constant is the Alfvén frequency, \omega_{\mathrm{A}}=k_{x}c_{\mathrm{A}}.
The black lines on the left panel of Figure 4 correspond to the predictions of the dispersion relation derived in Paper II for the linear Alfvénic modes and are the same as those represented in the top panel of Figure (2) of Paper II: the solid line corresponds to the Rmode and the dashed line represents the Lmode. We refer readers to Paper II for detailed explanations of the differences between these two solutions, which are only distinct from the classic Alfvén wave for high enough frequencies. The symbols represent the results of the numerical simulations. It can be seen that the simulations are in good agreement with the predictions of the dispersion relation for the firstorder perturbations.
Regarding the secondorder modes, their frequency, denoted by \omega_{sim}^{(2)}, has three clearly different regimes depending on the wavenumber of the firstorder perturbation. At small wavenumbers, the frequency is approximately given by \omega_{sim}^{(2)}\approx 2\widetilde{c}_{S}k_{x} (see Equation (A25)). However, as k_{x} is increased, \omega^{(2)} departs from that value. To understand this behavior, we take into account that at small wavenumbers, the oscillation frequency of the firstorder perturbation, \omega_{sim}^{(1)}, is lower than all collision frequencies. Thus, there is a considerably strong coupling between the three components of the plasma and they behave almost as a single fluid whose effective sound speed is given by \widetilde{c}_{S}. The resulting acoustic mode has a normalized frequency given by \omega_{sim}^{(2)}/\omega_{\mathrm{A}}=2\widetilde{c}_{S}k_{x}/\omega_{\mathrm% {A}}\approx 0.14. At intermediate wavenumbers, \omega_{sim}^{(1)} is larger than \nu_{p\mathrm{He}} and \nu_{\mathrm{HHe}}, but smaller than \nu_{p\mathrm{H}}, which means that neutral helium is weakly coupled to the other two fluids but protons and neutral hydrogen still have a strong interaction. Consequently, the effective sound speed is given by the weighted mean of those of protons and hydrogen, without the contribution of neutral helium, and is slightly larger than \widetilde{c}_{S}. With this new effective sound speed, the normalized oscillation frequency is \omega_{sim}^{(2)}/\omega_{\mathrm{A}}\approx 0.15. Finally, at large wavenumbers, \omega_{sim}^{(1)}\gg\nu_{p\mathrm{H}} and the coupling between protons and hydrogen is weak. The sound speed of the proton fluid is c_{ie} and \omega_{sim}^{(2)}/\omega_{\mathrm{A}}\approx 0.18, which corresponds approximately to the result expected for a fully ionized plasma. The neutral hydrogen and neutral helium fluids oscillate with the normalized frequencies 2k_{x}c_{S,\mathrm{H}}/\omega_{\mathrm{A}}\approx 0.13 and 2k_{x}c_{S,\mathrm{He}}\approx 0.065, respectively.
Furthermore, the simulations represented in Figure 13 reveal another contrast between the partially ionized and the fully ionized cases. In the latter, the ponderomotive force produces an accumulation of matter around the nodes of the Alfvén wave magnetic field perturbation. In a pressureless fluid, the accumulation continues without limit. However, when the effect of the gas pressure is taken into account, the density at that node reaches a certain maximum and starts to oscillate between that maximum and its background value (see, e.g., Rankin et al., 1994; Tikhonchuk et al., 1995). However, Figures 2 and 3 show that in partially ionized plasmas, the density tends to accumulate at the node only for a brief period of time. Later, the relative variation of density, \Delta\rho/\rho_{0}, decreases and oscillates around negative values, which means that the plasma becomes lighter at that location. This behavior can be understood in terms of the effect of gas pressure and collisions as follows.
The secondorder longitudinal motion of ions mainly depends on the balance between the forces given by the gradients of the thermodynamic and magnetic pressures. The study of fully ionized plasmas (see, e.g., Rankin et al., 1994; Tikhonchuk et al., 1995) shows that the gradient of magnetic pressure moves the plasma towards the nodes of the firstorder magnetic field perturbation. The gradient of the secondorder perturbation of pressure acts in the opposite way, i.e., it displaces the matter from those locations. Which effect dominates during the first steps of the temporal evolution depends on the time scales associated with them. Under the physical conditions used in this investigation, the Alfvén frequency is higher than the frequency of sound waves, meaning that magnetic pressure has a smaller time scale than the thermodynamic pressure. Therefore, initially, the matter accumulates at the nodes. Later, the effect of the thermodynamic pressure becomes noticeable and the resulting motion is a consequence of the combination of the two forces. In partially ionized plasmas, friction due to ionneutral collisions dissipates the energy of Alfvén waves and turns it into internal energy of the plasma, i.e., it increases the thermodynamic pressure. As time advances the term of the motion equation associated with the driving Alfvénic wave becomes less relevant in comparison with the gradient of the thermodynamic pressure. Consequently, the longitudinal motion is dominated by the force that moves the matter away from the nodes of the magnetic field perturbation.
Up to now the amplitude of the perturbations has been chosen in a way that only first and secondorder effects are relevant for the dynamics. However, if the amplitudes are increased, higherorder terms may be also of great importance. As detailed by Tikhonchuk et al. (1995), the higherorder terms may produce, for instance, the steepening of the fluctuations, which may lead to the formation of shocks, and the appearance of higher harmonics of the Alfvén waves. Some of those higherorder effects can be found in Figure 5, where the results of simulations with different amplitudes of the initial perturbation are compared. The steepening of the waves when the amplitudes are increased is clearly shown in panels (b) and (c), corresponding to the normalized xcomponent of the velocity at the point x=l/2 and the variation of density at x=0, respectively. The top left panel represents the firstorder Alfvén wave at x=0 and it can be seen that, after the initial steps, its frequency raises in the cases with the larger amplitudes. This is due to the decrease of density shown in panel (c). The change in frequency can also be noticed in the other three panels: a larger number of periods can be found for V_{y,0}=0.15c_{\mathrm{A}} than for V_{y,0}=0.025c_{\mathrm{A}}. Finally, panel (d) represents the variation of temperature at x=0. After a very fast growth, the temperature tends to oscillate around a value that increases with the square of the driver amplitude. This behavior is consistent with the heating term in Equation (1) of Paper II.
The results displayed in Figure 5(d) correspond to a specific point of the numerical domain. Although they are representative of the general heating of the plasma, differences appear (for example, in amplitude and in the phase of the oscillations) when other locations are considered. Thus, it is interesting to compute the average value over the spatial domain. The temporal evolution of the spatiallyaverage temperature, given by 1/(2l)\int_{l}^{l}T(x)\,dx, is represented in Figure 6. Comparing this figure with Figure 5, it can be seen that the temperature reaches an equilibrium value after most of the energy of the Alfvén wave has been dissipated, while the contribution of the secondorder acoustic waves to heating is negligible. This is due to ionneutral collisions being more efficient in damping the Alfvénic modes than in damping the acoustic modes under the parameters chosen for these simulations. The reason for this behavior is that collisional damping is more efficient when \omega\approx\nu_{st} (Leake et al., 2005; Zaqarashvili et al., 2011; Soler et al., 2013b) and in these simulations the frequency of the Alfvénic waves is closer to the collision frequencies than the frequency of the acoustic modes. When the amplitude of the initial perturbation is V_{y,0}=0.025c_{\mathrm{A}}, the temperature rises up to \sim 10,360\ \mathrm{K} (i.e., the variation is \Delta T\approx 360\ \mathrm{K}). For the amplitudes V_{y,0}=0.05c_{\mathrm{A}}, V_{y,0}=0.1c_{\mathrm{A}} and V_{y,0}=0.15c_{\mathrm{A}}, the final temperatures are 11,470\ \mathrm{K} (\Delta T\approx 1470\ \mathrm{K}), 16,170\ \mathrm{K} (\Delta T\approx 6170\ \mathrm{K}), and 24470\ \mathrm{K} (\Delta T\approx 14,470\ \mathrm{K}), respectively. Hence, the dependence of the increment of temperature on the amplitude of the perturbation is approximately quadratic.
3 Propagating waves: impulsive driver
In this section, the evolution of a velocity pulse as it propagates through a uniform partially ionized plasma is analyzed. A similar study was performed by Verwichte et al. (1999) for the case of fully ionized plasma. Hence, it is interesting to examine how the results of that work are modified by the inclusion of partial ionization effects. Moreover, according to Rankin et al. (1994), the effects of nonlinearity are stronger for standing waves than for propagating waves. Thus, third or higherorder terms are only expected to have a strong impact on the evolution of the pulse for larger amplitudes than those used in the previous section and the main nonlinearities that appear in this section are due to the secondorder terms. As before, we consider 1.5D numerical simulations and use physical conditions representative of solar prominence cores.
Now, the perturbation applied to the plasma at t=0 has a Gaussian profile, i.e., it is given by
f^{(1)}(x,t=0)\sim\exp\left[\left(\frac{xx_{0}}{\sqrt{2}\sigma_{x}}\right)^{% 2}\right],  (4) 
where \sigma_{x} is the rootmeansquare width and is related to the full width at half maximum (FWHM) of the Gaussian by the formula \mathrm{FWHM}=2\sqrt{2\ln 2}\sigma_{x}, and x_{0} is the central position of the peak.
Figure 7 shows the Alfvén wave that is generated when the perturbation given by Equation (4) is applied to the ycomponent of the velocity of all species. The amplitude of the perturbation is V_{y,0}=5\times 10^{2}c_{\mathrm{A}} and its width is \mathrm{FWHM}=2\times 10^{5}\ \mathrm{m}. As expected, the initial pulse splits into two smaller Alfvénic pulses, with half the height of the initial pulse, and propagate towards opposite directions. There is a strong coupling between the three species (protons, neutral hydrogen and neutral helium) and the transverse velocity pulses of each fluid propagate together at the modified Alfvén speed, \widetilde{c}_{\mathrm{A}}. Notwithstanding, the height of the peaks decreases with time because the coupling is not perfect and there is friction that dissipates a fraction of the wave energy and turns it into internal energy of the plasma. Friction is caused by the small velocity drifts between species, which are not noticeable at the scale of Figure 7.
The nonlinear effects generated by the Alfvénic pulse are represented in Figure 8. The panels in the top row display the perturbation on the xcomponent of the velocity. The amplitude of V_{x} is much smaller than that of V_{y}, of the order of 1.5\% of V_{y,0}, as it would be expected. As in the case of standing waves, two clearly different propagating waves appear in the longitudinal component of velocity. The faster one has a propagation speed that coincides with \widetilde{c}_{\mathrm{A}}, while the slower one propagates at the speed \widetilde{c}_{S}. The waves leave a small wake that is positive at x>0 and negative at x<0. This means that, after the wavefront has passed, the particles are slowly moved away from the center. Again, this is a nonlinear effect.
The relative variation of density is shown in the second row of Figure 8. Although their shapes are different, the perturbations found here have the same propagation speeds as those for V_{x}. Moreover, a similar behavior to that previously described for standing waves can be observed: matter accumulates at the center of the domain during the first steps of the simulation but later is displaced from that location.
The third and fourth rows of Figure 8 represent the secondorder perturbations of pressure and temperature, respectively, with \Delta P=P(x,t)P_{0}. These two rows show how a fraction of the energy of the perturbation is deposited into the plasma. An increase of temperature and pressure is found after the passing of the wave front, i.e., some of the energy of the wave has been transformed into internal energy of the plasma. The increase of pressure seems to be uniform along the plasma. In contrast, it can be checked that the growth of temperature is inversely proportional to the variation of density.
The results shown in Figure 7 and the first and second columns of Figure 8 can be compared with those in Figure 1 from Verwichte et al. (1999). A similar behavior is found in fully and partially ionized plasmas during the first steps of the evolution of the density and velocity. The differences appear in pressure and in temperature. Verwichte et al. (1999) did not plot the evolution of the pressure because, for the case of fully ionized plasmas, it has the same shape as that of density. In contrast, in partially ionized plasmas, the propagating waves leave a pressure wake due to the frictional dissipation of energy because of ionneutral collisions, a phenomenon that is obviously absent from the fully ionized case of Verwichte et al. (1999).
Not all the kinetic energy of the initial perturbation is used in heating the plasma, but a fraction of it is inverted in generating the secondorder propagating waves. Hence, it is interesting to investigate how the energy deposition depends on the properties of the initial perturbation. A series of simulations have been performed with different widths of the Gaussian velocity pulse but keeping the same initial kinetic energy, i.e., the amplitude of the pulse has been modified accordingly. The results of this study are displayed in Figure 9.
The background internal energy is computed after the two wavefronts, i.e., the Alfvénic pulse and the nonlinearly generated sonic pulse, have abandoned the the numerical domain of interest. The initial kinetic energy is computed as
e_{k}(t=0)=\frac{1}{2l}\int_{l}^{l}\sum_{s}\rho_{s}(x,t=0)\left[V_{y}^{(1)}(x% ,t=0)\right]^{2}\,dx,  (5) 
and the variation of the internal energy of the medium is given by
\Delta e_{P}(t)=\frac{1}{2l}\int_{l}^{l}\sum_{s}\frac{P_{s}^{(2)}(x,t)P_{s,0% }(x)}{\gamma1}\,dx.  (6) 
Left panel of Figure 9 shows that the deposition of energy into the plasma has a remarkable dependence on the width of the pulse. A peak of \Delta e_{p}/e_{k}(t=0)\approx 6\% is found at \mathrm{FWHM}=10^{5}\ \mathrm{m}, which corresponds to a perturbation with an amplitude of V_{y,0}=0.1c_{\mathrm{A}}/\sqrt{2}. At larger widths, the fraction of deposited energy decreases exponentially. This behavior can be understood by taking into account that the width of a Gaussian pulse is associated with a certain scale of wavelengths or wavenumbers. Perturbations with larger widths are associated to smaller scales of wavenumbers and at smaller wavenumbers the coupling between the species of the plasma is stronger and the dissipation of energy is smaller.
Additional series of simulations have been performed to check if the trend examined in the previous paragraphs is also found under different conditions. In the first set of new simulations we apply the initial perturbation only to the ions, leaving neutrals initially at rest. In another series of simulations we perturb the ycomponent of the magnetic field instead of the velocity. The results are represented in the middle and right panels of Figure 9, respectively. For the latter case, the magnetic energy density of the initial perturbation has been computed as
e_{B}(t=0)=\frac{1}{2l_{1}}\int_{l_{1}}^{l_{1}}\frac{\left[B_{y}^{(1)}(x,t=0)% \right]^{2}}{2\mu_{0}}\,dx.  (7) 
The comparison of the left and middle panels of Figure 9 shows the same type of dependence of the energy deposition on the width of the perturbation. However, the peak value is \sim 2\% when neutrals are initially at rest instead of \sim 6\% when the perturbation is applied to all species. The reason may be that a considerable fraction of the energy has to be used in setting the neutrals in motion by means of collisions with ions: it must be reminded that under the chosen prominence conditions, neutrals account for 2/3 of the total mass of the plasma.
When the perturbation is applied to the ycomponent of the magnetic field (right panel of Figure 9), the dependence of the energy transfer is similar to the one found in the previous cases. The peak appears at \mathrm{FWHM}\approx 10^{5}\ \mathrm{m} and it has the same value as in the left panel of Figure 9, \sim 6\%. So, regarding the eventual energy deposition into the plasma due to wave dissipation, it is irrelevant whether the energy of the initial perturbation is kinetic or magnetic, as long as the total energy is the same.
The results described in the paragraphs above seem to be in good agreement with the findings of Papers I and II, i.e., larger wavenumbers are more damped than smaller ones. However, for very small values of the perturbation width, Figure 9 shows a peculiar trend that diverges from what it might be expected: the efficiency of energy deposition decreases as the width of the initial perturbation is reduced (and the associated wavenumbers are larger). The reason may be related to the the fact that quite large amplitudes of the perturbations are needed when the widths are reduced in order to keep the initial energy the same in all simulations. As already mentioned, the energy of the initial perturbation is used in two ways, namely generation of waves and heating of the plasma. Hence, the internal energy has two components: one associated to the propagating wavefronts and another one related to energy gains and losses of the background plasma. A study of how those two components vary is illustrated by Figure 10, where the temporal evolution of the kinetic, magnetic, internal and total energy is displayed for four simulations.
In Figure 10 the total energy is not constant but diminishes with time: the waves are leaving the region of interest, carrying with them an important fraction of the initial total energy. This can be clearly noticed at t\approx 10\ \mathrm{s}, when most of the kinetic and magnetic energy goes to zero because Alfvén waves start crossing the boundaries. Later, the nonlinearlygenerated sound waves also abandon the domain and the remaining energy is, then, truly associated with what is deposited in the plasma.
It must be noted that the peak that can be seen at t\approx 80\ \mathrm{s} is a consequence of the sound waves leaving the domain of interest. It does not mean that there is a sudden increase of energy in the simulation: the total energy remains constant if we also account for the energy of the escaped waves. The peak appears because the leading section of the sound wave has a negative contribution to the perturbation of the internal energy (as can be seen in the third row of Figure 8) and, as it leaves the domain, generates the effect of an apparent rise of energy.
Focusing on the first seconds of the simulations, it can be seen that the amount of the initial energy that is transformed into internal energy increases with the amplitude of the perturbation (or, equivalently, when the width diminishes): the height of the dashed line (which represents the internal energy) at t\approx 10\ \mathrm{s} is larger in the bottom right panel, which corresponds to an amplitude of V_{y,0}=0.2c_{\mathrm{A}} and \mathrm{FWHM}=1.25\times 10^{4}\ \mathrm{m}. Thus, a larger amplitude of the initial perturbation corresponds to a larger increment of the internal energy. However, the distribution of this increment between the energy associated to the propagating wavefronts and that actually deposited into the plasma is not always the same: for instance, although the increase of internal energy is larger for V_{y,0}=0.2c_{\mathrm{A}} than for V_{y,0}=0.1c_{\mathrm{A}} (bottom left panel), at the end of the simulation the latter case retains more internal energy. This means that the contribution from waves represents a larger fraction of the internal energy when the amplitude of the perturbation increases, i.e., when the nonlinear effects are more relevant. The reason is that more energy is required to generate secondorder waves when the amplitude of the firstorder perturbation increases, which leaves a smaller fraction of the initial energy that can be used in heating the plasma.
4 Propagating waves: periodic driver
After analyzing the properties of nonlinear Alfvén waves generated by an impulsive driver, we turn to the case of waves excited by a periodic driver. Here, we perform numerical simulations in which a linearly polarized harmonic driver is applied to the boundary x=0 and the resulting waves are allowed to propagate along the positive xaxis. We again consider physical conditions that correspond to a quiescent prominence and the driver is given by a periodic function of time. As a boundary condition, we impose that V_{x,s}(x=0,t)=0 to avoid any mass inflow through the boundary from outside the numerical domain. All other variables are extrapolated at x=0.
Figure 11 shows the results of a simulation where the driver is applied to the ycomponent of the velocity of all species during a time t_{D}, after which the driver is switched off. Hence, the driver can be cast as
V_{y,s}(x=0,t)=\left\{\begin{array}[]{l c}V_{y,0}\sin\left(\omega t\right)&% \text{if $t\leq t_{D}$},\\ 0&\text{if $t>t_{D}$}.\end{array}\right.  (8) 
where the amplitude of the driver is V_{y,0}=0.1c_{\mathrm{A}} and its frequency is \omega=10^{5}\Omega_{p}\approx 1\ \mathrm{rad\ s^{1}}, with \Omega_{p}=Z_{p}eB_{0}/m_{p} the cyclotron frequency of protons and Z_{p} the charge number. In this case, the driver has been applied during 3 periods of the Alfvén wave, i.e., t_{D}=3\tau, where \tau=2\pi/\omega. In the same way as in the case of standing waves, the driver used here excites circularly polarized waves, which means that perturbations in the other transverse component, i.e., the zcomponent, of the velocity and the magnetic field are also generated. However, due to the small amplitudes of these perturbations, in this section we again focus only on the ycomponents.
The figure displays the results after the driver has been switched off. The top left panel shows the driven firstorder wave, which propagates at the modified Alfvén speed, \widetilde{c}_{\mathrm{A}}. The dotted lines in this panel correspond to the spatial damping predicted by Equation (16) of Paper II for linear perturbations. The numerical results are in good agreement with that prediction.
The top right panel of Figure 11 shows the nonlinearly generated perturbations in the longitudinal component of the velocity. Two waves can be clearly identified. The fastest one travels with velocity \widetilde{c}_{\mathrm{A}} and its longitudinal velocity is always positive. The slowest perturbation is associated with the effective sound speed, \widetilde{c}_{S}, and its longitudinal velocity is always negative. While the driver is on, the longitudinal motion of the plasma in the region x\leq\widetilde{c}_{S}t results from the combination of the two waves. The frequency of the two waves is twice the frequency of the firstorder Alfvén wave and their wavenumbers are also larger. A similar behavior can be seen in the bottom left panel, where the relative variation of density is plotted. The density decreases in the region x\leq\widetilde{c}_{S}t while it increases in \widetilde{c}_{S}<x<\widetilde{c}_{\mathrm{A}}t. This means that the driver creates a mass flow that displaces the plasma away from the boundary. The aforementioned characteristics of the secondorder perturbations are consistent with the analytical results shown in Appendix A.2, which have been derived for the simpler case of a twofluid partially ionized plasma with strong ionneutral coupling.
The bottom right panel of Figure 11 shows the relative variation of temperature. This variation comes from a combination of two effects: the fluctuations of density and pressure associated with the nonlinearly generated longitudinal perturbations and the rise of the internal energy of the plasma due to the collisional friction. Although the driver used in this simulation fulfills \omega/(2\pi)<\nu_{st} and the coupling between the components of the plasma is strong, there are still small velocity drifts, not noticeable at the scale of the figure, which lead to the dissipation of a fraction of the energy of the driver. Consequently, the temperature of the plasma increases as the perturbations propagate. The largest growth of temperature occurs near the boundary where the driver is applied. The temperature gradients create, in turn, pressure gradients that further contribute to displace the plasma from the boundary. The effect of this pressure gradient can be clearly seen, for instance, on the left side of the region \widetilde{c}_{S}t<x<\widetilde{c}_{\mathrm{A}}t of the top right panel, where V_{x,s}(x,t)>0, indicating that there is a mass flow towards the right.
Next, we study how the amplitude of the driver affects the properties of the nonlinear generated waves and the growth of the internal energy of the plasma. From what has been found in the previous sections, it would be expected that the total energy deposited in the plasma has a quadratic dependence on the amplitude of the driver. The top left panel of Figure 12 shows the relative variation of internal energy in three simulations with V_{y,0}=0.025c_{\mathrm{A}}, V_{y,0}=0.05c_{\mathrm{A}}, and V_{y,0}=0.1c_{\mathrm{A}}, respectively. There is a slight decrease of the internal energy around t=40\ \mathrm{s}: at that time, the Alfvénic perturbations start to leave the domain. It can be checked that the relative variation of internal energy for the case with V_{y,0}=0.1c_{\mathrm{A}} is approximately four times larger than the relative variation for the case with V_{y,0}=0.05c_{\mathrm{A}}, and approximately sixteen times larger that for V_{y,0}=0.025c_{\mathrm{A}}, which corresponds to a quadratic dependence on the amplitude of the driver, as expected.
The right panel shows the perturbations of the longitudinal component of the velocity of ions at a given time of the simulations. The plot focus on a region close to the boundary to better examine the shape of the perturbations. Increasing the amplitude of the driver produces a steepening of the waves, which is more pronounced for sound waves than for Alfvénic waves. The reason is that, under the chosen physical conditions, the effective sound speed is smaller than the modified Alfvén speed and, hence, the amplitude of the driver is highly nonlinear if compared with \widetilde{c}_{S} but it is not so large if compared with \widetilde{c}_{\mathrm{A}}. Consequently, sound waves develop shocks more easily than Alfvénic waves.
The development of shocks is related to the variation of the propagation speeds with respect to their values at the equilibrium state. This variation is represented in Figure 13, where only the results of the simulations with V_{y,0}=0.1c_{\mathrm{A}} and V_{y,0}=0.05c_{\mathrm{A}} are shown. As the leading perturbation travels to the right, the propagation speeds of the trailing waves change. One of the reasons of this change is the fluctuations in density. For instance, in the region \widetilde{c}_{S}(t=0)t<x<\widetilde{c}_{\mathrm{A}}(t=0)t the relative variation of density is positive, as shown in Figure 11, and the Alfvén speed is smaller than in the region x>\widetilde{c}_{\mathrm{A}}(t=0)t. On the other hand, \Delta\rho/\rho_{0}<0 for x<\widetilde{c}_{S}(t=0)t and, consequently, \widetilde{c}_{\mathrm{A}}(x,t)>\widetilde{c}_{\mathrm{A}}(t=0). The propagation speed of the sound waves is also modified by those fluctuations on density but it is also affected by the variation of pressure. A larger pressure implies a larger effective sound speed. In a partially ionized plasma, the increase of the effective sound speed is enhanced by the dissipation of the energy of the firstorder Alfvén wave due to ionneutral collisions. Therefore, the secondorder acoustic waves turn into shocks more easily in partially ionized plasmas than in fullyionized plasmas. Figure 13 shows that the fluctuations of \widetilde{c}_{\mathrm{A}} are smaller than those of \widetilde{c}_{S}, which explains the differences in the steepening of the Alfvénic and the sound waves represented in the right panel of Figure 12.
The formation of shocks through the ponderomotive coupling of Alfvén waves to sound modes was investigated by Arber et al. (2016). Their 1.5D numerical study suggests that in the chromosphere the heating due to shocks is larger than that caused directly by ionneutral collisions. Here, we have shown that in a quiescent prominence sound waves develop shocks more easily than Alfvén waves. Hence, shock heating may have some contribution to the total heating of partially ionized prominences. Nevertheless, since viscosity is not included in our model, we cannot compute the heating associated with the dissipation of acoustic shocks.
Now, we perform simulations in which the driver is active during all the running time of the simulation and not only for a brief time interval. Our goal here is to study the heating rate due to ionneutral collisions. For this series of simulations the driver is given by
B_{y,s}(x=0,t)=B_{y,0}\sin\left(\omega t\right),  (9) 
with B_{y,0}=10^{2}B_{0}. Thus, the energy input by the effect of the driver is given by
e_{B}=\frac{B_{y}(x=0,t)^{2}}{2\mu_{0}}=\frac{B_{y,0}^{2}}{2\mu_{0}}\sin^{2}% \left(\omega t\right)  (10) 
and the mean energy input per period is
\langle e_{B}\rangle=\frac{1}{\tau}\int_{0}^{\tau}e_{B}\,dt=\frac{1}{4}\frac{B% _{y,0}^{2}}{2\mu_{0}},  (11) 
Note that the mean energy is independent of the frequency of the driver. The reason for exciting the waves through the perturbation of the magnetic field instead of the velocity as before is that we require the energy input to be constant. Due to the density variations caused by the nonlinear waves, the energy input is not constant if the driver is applied to the velocity as before.
From the mean energy density, we compute the power of the driver as
\mathcal{P}=\frac{\langle e_{B}\rangle}{\tau}=\frac{\omega\langle e_{B}\rangle% }{2\pi},  (12) 
which has the same physical units as the heating rate, namely W\ m^{3}, and, hence, the two quantities can be directly compared. The formula above shows that for a fixed amplitude of the driver, perturbations of higher frequency have a larger power.
The heating rate can be computed as
Q=\sum_{s}\sum_{t\neq s}Q_{s}^{st},  (13) 
where Q_{s}^{st} is given by Equation (1) of Paper II. For the case of the threefluid prominence plasma, the heating rate is then given by
\displaystyle Q  \displaystyle=  \displaystyle\alpha_{p\mathrm{H}}\left(\bm{V_{p}}\bm{V_{\mathrm{H}}}\right)^{% 2}+\alpha_{p\mathrm{He}}\left(\bm{V_{p}}\bm{V_{\mathrm{He}}}\right)^{2}+% \alpha_{\mathrm{HHe}}\left(\bm{V_{\mathrm{H}}}\bm{V_{\mathrm{He}}}\right)^{2}  (14)  
\displaystyle+  \displaystyle\frac{\alpha_{pe}}{(en_{e})^{2}}\bm{j}^{2}+\alpha_{e\mathrm{H}}% \left(\bm{V_{\mathrm{H}}}\bm{V_{p}}+\frac{\bm{j}}{en_{e}}\right)^{2}+\alpha_{% e\mathrm{He}}\left(\bm{V_{\mathrm{He}}}\bm{V_{p}}+\frac{\bm{j}}{en_{e}}\right% )^{2}, 
where \bm{j}=\sum_{s}Z_{s}en_{s}\bm{V_{s}} is the current density and the last three terms correspond to the effect of collisions with electrons, i.e., magnetic resistivity.
In the first place, we examine the spatial distribution of the heating rate. Figure 14 shows the results of a simulation with a driver of amplitude B_{y,0}=10^{2}B_{0} and frequency \omega=2\times 10^{4}\Omega_{p}. The vertical dotted lines represent the distance at which the amplitude of the wave has been reduced by a factor \exp(3) due to collisional damping. This distance is given by \lambda_{I}=3/k_{I}, where k_{I} is the imaginary part of the solution to the dispersion relation for linear Alfvén waves given by Equation (16) from Paper II. The choice of this reference distance is somewhat arbitrary but it will be useful in forthcoming computations.
The top panel of Figure 14 shows the normalized ycomponent of the magnetic field. The numerical results are in good agreement with the damping of the wave predicted by the dispersion relation: as the wave propagates along the xaxis, its amplitude is proportional to \exp(k_{I}x).
The bottom panel shows the heating rate. Its amplitude decays with the distance much faster than the magnetic field. Fitting the heating rate with an exponentially decaying function, we check that it is proportional to \exp(2k_{I}x). This is the expected behavior since, according to Equation (14), Q has a quadratic dependence.
Next, we study the dependence of the heating rate on the frequency of the driver. Thus, we keep the amplitude of the driver fixed at B_{y,0}=10^{2}B_{0} and perform a series of simulations with different frequencies. Then, we compute the spatiallyaveraged heating rate at every time step of the simulation in the following way:
Q_{av}(t)=\frac{1}{\lambda_{I}}\int_{0}^{\lambda_{I}}Q(x,t)\,dx.  (15) 
During the first time steps, Q_{av} grows as the waves propagate along the xaxis. After the waves reach the position x=\lambda_{I}, the profile of Q_{av}(t) flattens and it tends to a constant value with small fluctuations. Hence, we compute Q_{av} as the mean value of Q_{av}(t) without taking into account the initial time steps. The results of this study are shown in Figure 15.
Song & Vasyliūnas (2011) found that the heating has a quadratic dependence on \omega when \omega\ll\nu_{st} and that it is independent from \omega if the oscillation frequency is much larger than the collision frequencies. On the contrary, Figure 15 suggests three different dependencies. To check that, the results are fitted to the powerlaw function Q_{fit}\sim\omega^{a}. In the range of low frequencies, the fitting exponent is a\approx 1.93, quite close to the quadratic dependence found by Song & Vasyliūnas (2011). The exponent for intermediate frequencies switches to a\approx 1.32. Finally, for higher oscillation frequencies, heating becomes independent of \omega. The differences between our results and those of Song & Vasyliūnas (2011) reside in the number of considered species. Song & Vasyliūnas (2011) only considered one ionized and one neutral species. Hence, at very low frequencies, \omega is much lower than all the collision frequencies. Therefore, all the species are strongly coupled and heating tends to a quadratic dependence on \omega, as Song & Vasyliūnas (2011) found. However, as \omega increases, the resulting behavior is a combination of the different trends associated with the various collision frequencies and departs from the quadratic dependence of the simplified twofluid case of Song & Vasyliūnas (2011). Finally, at very high oscillation frequencies, heating does not depend on \omega, in agreement with Song & Vasyliūnas (2011).
Figure 15 also allows to compare the heating rate with the power of the driving waves. This gives a measure of the efficiency of collisions as a heating mechanism. At low oscillation frequencies, the heating rate is much lower than the input power. This means that only a very small fraction of the energy added to the plasma by the driver is used in increasing the temperature of the fluid. However, this fraction increases at larger frequencies but it decreases again when heating becomes independent of the oscillation frequency. Consequently, heating by ionneutral collisions is more efficient when \omega\approx\nu_{st} (Leake et al., 2005; Zaqarashvili et al., 2011; Soler et al., 2013b).
If the driver acts continuously on the plasma, the temperature rises without bound. However, in reality radiative losses of energy would prevent this unlimited increase. It is interesting to check if the heating rates computed are large enough to balance radiative losses, which will be denoted by L. From Parenti et al. (2006), Parenti & Vial (2007), and Heinzel et al. (2010), Soler et al. (2016) estimated that the volumetric radiative losses of a prominence are of the order of 10^{5} to 10^{4}\ \mathrm{W\ m^{3}}. Comparing those values with the results represented in Figure 15, we find that L<Q for \omega<10^{4}\ \Omega_{p}. At low frequencies, heating by ionneutral collisions can only compensate a fraction of the radiated energy. In the range between \omega\approx 10^{4}\ \Omega_{p} and \omega\approx 10^{3}\ \Omega_{p}, Q becomes of the order of L. At frequencies higher than \omega\approx 10^{3}\ \Omega_{p}, the heating rates are larger than the radiative losses. Nevertheless, it must be taken into account that at high frequencies the energy of the waves is dissipated in a smaller region than at lower frequencies. Hence, collisions would only heat that region of the prominence. In addition, we remind that the results of Figure 15 have been obtained for a given driver amplitude and heating rates strongly depend on the amplitude of the driving wave.
5 Conclusions
In this paper, we have continued the investigation started in Papers I and II about the importance of multifluid effects for the correct description of waves in multicomponent, partially ionized plasmas. In Papers I and II we focused on the study of linear waves. In the present and final installment of the series, we have addressed nonlinear effects, paying special attention to the role of particle collisions in plasma heating.
Nonlinear waves in partially ionized plasmas have been studied in this work by means of a multifluid model in which the effects of elastic collisions between all species of the plasma are taken into account. The general properties of nonlinear lowfrequency Alfvén waves analyzed here are consistent with the results obtained by, e.g., Hollweg (1971), Rankin et al. (1994), Tikhonchuk et al. (1995) or Verwichte et al. (1999) for fully ionized plasmas, although differences appear due to the collisional interaction between ions and neutrals. For example, a secondorder effect of nonlinear standing Alfvén waves is the appearance of a ponderomotive force that induces fluctuations in density, pressure and the longitudinal component of the velocity. For the case of standing waves in fully ionized plasmas, those variations are a combination of two modes with frequencies given by 2k_{z}c_{\mathrm{A}} and 2k_{z}c_{\mathrm{ie}}, and their wavenumber is twice the value for the original perturbation. However, in partially ionized plasmas, the frequencies are proportional to the modified Alfvén speed, \widetilde{c}_{\mathrm{A}} and the weighted mean sound speed, \widetilde{c}_{S}, respectively, when the smallwavenumber range is considered. Since in the plasmas that have been examined here \widetilde{c}_{\mathrm{A}} is much lower than \widetilde{c}_{S}, the secondorder oscillations induced by the ponderomotive force are dominated by the mode associated with the sound speed. Due to this ponderomotive force, the matter of the plasma tends to accumulate at the nodes of the magnetic field wave, although such accumulation is limited by the effect of pressure.
If the wavenumber of the perturbations is increased, the coupling between the different species is reduced and the collisional friction becomes relevant. It is then that multifluid effects become of interest. The plasma is heated and the effect of pressure against the accumulation of matter is enhanced. After the original Alfvén wave has dissipated due to collisions, the result of the action of the ponderomotive force is the displacement of matter from the nodes of the magnetic field towards the antinodes instead of its accumulation at the nodes. At even higher frequencies, the species of the plasma become almost uncoupled from each other and the oscillation frequencies of the secondorder waves tend to the values predicted for fully ionized plasmas, although they are strongly damped because of collisions. These results were obtained through the study of an initial perturbation that was weakly nonlinear. Cases with larger amplitudes have also been briefly analyzed and it was found that the profile of the nonlinear waves steepens as time advances and the frequency of the oscillations are slightly modified due to the more important variations of density, which is consistent with the findings of Tikhonchuk et al. (1995) or Verwichte et al. (1999).
The propagation of nonlinear pulses through a plasma with conditions akin to those of a solar quiescent prominence has also been examined. The simulations have shown that after the initial perturbation has been applied to the plasma, the pulse splits in two smaller pulses that propagate in opposite directions at a speed given by \widetilde{c}_{\mathrm{A}}. The amplitude of those pulses decreases with time due to the collisions between ions and neutrals, which dissipates a fraction of the energy of the initial perturbation. The amount of dissipated energy increases when the width of the perturbation decreases. This behavior is due to the larger wavenumbers associated with a smaller width of the pulse. According to the results from Paper II, waves with larger wavenumbers have shorter damping times due to ionneutral collisions, while perturbations with smaller wavenumbers are more longlived. Hence, the widths of the pulses increase and their amplitudes diminish with time as the larger wavenumbers are dissipated by the collisional friction.
As a secondorder effect, the pulse generates two pairs of longitudinal waves that propagate in opposite directions. The phase speeds of those waves are given by \widetilde{c}_{\mathrm{A}} and \widetilde{c}_{S}, respectively, as one of them is associated to the primary Alfvén wave and the other one with a nonlinearly generated sound wave. In addition, a fraction of the initial energy is deposited in the plasma in form of heat. Consequently, the temperature of the plasma rises. The numerical simulations show that the heating generally increases when the width of the pulse is decreased, which, as already mentioned, is associated with the efficient dissipation of small scales. However, at small enough widths, the computed heating decreases again. This may be explained by the highly nonlinear amplitude of the perturbations. When the amplitude of the initial perturbation is increased, the generation of the secondorder waves requires a larger fraction of the initial energy; hence, there is a smaller fraction of energy available to be transformed into heat from the firstorder wave. For conditions of solar quiescent prominences, the investigation presented here has found that a maximum of a 6\% of the energy of the initial perturbation is finally used in heating the plasma. However, this value may vary if longer times, larger domains or different physical conditions are chosen for the simulations.
We have also studied the properties of nonlinear Alfvén waves generated by a periodic driver. As in the case of the propagating pulse, the driven Alfvén wave generates two secondorder perturbations in density, pressure and the longitudinal component of the velocity, which propagate at the speeds \widetilde{c}_{S} and \widetilde{c}_{\mathrm{A}}. As time advances, those perturbations cause a decrease of the density in the region x\leq\widetilde{c}_{S}t and an increase in the region \widetilde{c}_{S}t<x\leq\widetilde{c}_{\mathrm{A}}t. This means that the plasma is displaced towards the direction of propagation of the nonlinear Alfvén wave. In a partially ionized plasma, the gradient of pressure caused by the collisional friction also contributes to this effect. In addition, we have shown that ionneutral collisions enhance the formation of shocks from the secondorder sound waves.
We have computed the heating rates due to elastic collisions and studied its spatial distribution and its dependence on the frequency. On the one hand, if the amplitude of the driving wave decays as \exp(k_{I}x), the heating rate decays as \exp(2k_{I}x). On the other hand, when the oscillation frequency of the wave is much lower than the collisions frequencies, the heating rates tend to a quadratic dependence on \omega. In the opposite limit, the heating rates are independent of the oscillation frequency. A more complex dependence is obtained at the intermediate range of frequencies. At this range, some of the species of the plasma are weakly coupled to the rest, while others are still strongly coupled. Consequently, the dependence of heating on frequency is an intermediate state between both limits.
The heating rates have been compared with an estimation of the radiative losses in a prominence (L\sim 10^{5}10^{4}\ \mathrm{W\ m^{3}}). We have found that at low frequencies, collisional heating can only balance a small fraction of the radiative losses. Nevertheless, at higher oscillation frequencies, the heating rates are of the order of or larger than the radiative losses.
In the present work, we have focused on media that are initially homogeneous and static. The effects of inhomogeneities, such as the gravitationally stratified plasma of the solar chromosphere, and realistic geometries for solar prominences should be investigated in the future. Furthermore, here we have limited ourselves to the simplest case of 1.5D simulations. More realistic studies in 2D and 3D, which would allow us to explore in depth the properties of magnetoacoustic waves, are left for future works. This series of papers was meant to pave the way for more elaborated future studies that should exploit the full applicability of the nonlinear multifluid code developed here.
Appendix A Approximate analysis of nonlinear waves in a partially ionized twofluid plasma
Here, a partially ionized twofluid plasma is considered as a simpler, toy model that can help us to understand the numerical results given in the paper. One of the fluids is composed of ions and electrons, and the other one is composed of neutrals. For the sake of simplicity, Hall’s term and Ohm’s diffusion are neglected from the induction equation. Therefore, the equations that describe the dynamics of this plasma are a simplified version of those used in the full simulations, namely
\frac{\partial\rho_{i}}{\partial t}+\nabla\cdot\left(\rho_{i}\bm{V_{i}}\right)% =0,  (A1) 
\frac{\partial\rho_{n}}{\partial t}+\nabla\cdot\left(\rho_{n}\bm{V_{n}}\right)% =0,  (A2) 
\frac{\partial\left(\rho_{i}\bm{V_{i}}\right)}{\partial t}+\nabla\cdot\left(% \rho_{i}\bm{V_{i}}\bm{V_{i}}\right)=\nabla P_{ie}+\frac{\nabla\times\bm{B}}{% \mu_{0}}\times\bm{B}+\alpha_{in}\left(\bm{V_{n}}\bm{V_{i}}\right),  (A3) 
\frac{\partial\left(\rho_{n}\bm{V_{n}}\right)}{\partial t}+\nabla\cdot\left(% \rho_{n}\bm{V_{n}}\bm{V_{n}}\right)=\nabla P_{n}+\alpha_{in}\left(\bm{V_{i}}% \bm{V_{n}}\right)  (A4) 
and
\frac{\partial\bm{B}}{\partial t}=\nabla\times\left(\bm{V_{i}}\times\bm{B}% \right),  (A5) 
where P_{n} is the pressure of neutrals, P_{ie} is the sum of the pressures of ions and electrons and \alpha_{in} is the ionneutral friction coefficient. The expression of the friction coefficient can be found in Paper II. The rest of symbols have been defined before.
To study the properties of nonlinear waves, a perturbative expansion is performed. Thus, each variable, \bm{f}, in the previous system of equations is rewritten as follows:
\bm{f}=\bm{f}^{(0)}+\epsilon\bm{f}^{(1)}+\epsilon^{2}\bm{f}^{(2)}+\dots,  (A6) 
where \epsilon is a dimensionless parameter proportional to the velocity amplitude of Alfvén waves, the superscript “(0)” refers to the background values and the superscripts “(1)” and “(2)” correspond to the firstorder and secondorder perturbations, respectively. Since a static uniform background is considered, \bm{V_{i}}^{(0)}=\bm{V_{n}}^{(0)}=0 and the remaining background values are constant.
Then, the terms in Equations (A1)(A5) can be gathered according to their powers of \epsilon, and separated systems of equations can be obtained for each order of the perturbative expansion.
If the initial perturbations are chosen to be transverse to the direction of the background magnetic field (assumed here to be in the xdirection) and let them to propagate along that same direction, the firstorder (or linear) system leads to the equation for Alfvén waves,
\left[\frac{\partial^{3}}{\partial t^{3}}+\left(1+\chi\right)\nu_{ni}\frac{% \partial^{2}}{\partial t^{2}}c_{\mathrm{A}}^{2}\frac{\partial}{\partial t}% \frac{\partial^{2}}{\partial x^{2}}c_{\mathrm{A}}^{2}\nu_{ni}\frac{\partial^{% 2}}{\partial x^{2}}\right]\bm{V_{i,\bot}}^{(1)}=0,  (A7) 
where \chi=\rho_{n}/\rho_{i} is the ionization ratio, \nu_{ni}=\alpha_{in}/\rho_{n} is the neutralion collision frequency, and \bm{V_{i,\bot}}^{(1)}\equiv V_{i,y}^{(1)}\bm{\hat{\jmath}}+V_{i,z}^{(1)}\bm{% \hat{k}} is the perturbation of the velocity of ions in the perpendicular direction to the background magnetic field. The firstorder perturbation of magnetic field can be found through the equation
\frac{\partial\bm{B_{\bot}}^{(1)}}{\partial t}=B_{0}\frac{\partial\bm{V_{i,% \bot}}^{(1)}}{\partial x},  (A8) 
where B_{0}\equiv B_{x}^{(0)} is the background magnetic field and \bm{B_{\bot}}^{(1)}\equiv B_{y}^{(1)}\bm{\hat{\jmath}}+B_{z}^{(1)}\bm{\hat{k}}.
The solutions of Equation (A7) in the form of Fourier modes have been analyzed by, e.g., Piddington (1956), Kulsrud & Pearce (1969), Pudritz (1990), Martin et al. (1997), Kamaya & Nishi (1998), Kumar & Roberts (2003), Zaqarashvili et al. (2011), Mouschovias et al. (2011) or Soler et al. (2013b). At first order, there is no coupling between the perpendicular and longitudinal components of the perturbations, which means that there is no coupling between Alfvén and sound waves. In contrast, a coupling appears at the secondorder, as shown by the following equations, which are related to the velocities in the longitudinal direction:
\frac{\partial\rho_{n}^{(2)}}{\partial t}+\rho_{n}^{(0)}\frac{\partial V_{n,x}% ^{(2)}}{\partial x}=0,  (A9) 
\frac{\partial\rho_{i}^{(2)}}{\partial t}+\rho_{i}^{(0)}\frac{\partial V_{i,x}% ^{(2)}}{\partial x}=0,  (A10) 
\rho_{n}^{(0)}\frac{\partial V_{n,x}^{(2)}}{\partial t}=\frac{\partial P_{n}^% {(2)}}{\partial x}+\alpha_{in}\left(V_{i,x}^{(2)}V_{n,x}^{(2)}\right),  (A11) 
\rho_{i}^{(0)}\frac{\partial V_{i,x}^{(2)}}{\partial t}=\frac{\partial P_{ie}% ^{(2)}}{\partial x}\frac{\partial}{\partial x}\left(\frac{B_{\bot}^{2}}{2\mu_% {0}}\right)+\alpha_{in}\left(V_{n,x}^{(2)}V_{i,x}^{(2)}\right),  (A12) 
where B_{\bot}^{2}\equiv\left(B_{y}^{(1)}\right)^{2}+\left(B_{z}^{(1)}\right)^{2}. Thus, the secondorder perturbation of the velocity of ions is related to the firstorder perturbation of the magnetic field and, in turn, produces a fluctuation in the rest of the variables, namely V_{n,x}^{(2)}, \rho_{i}^{(2)}, and \rho_{n}^{(2)}. It must be noted that the secondorder equations corresponding to the perpendicular components have the same form as those of firstorder and hence, they describe the same behavior as Equations (A7) and (A8).
The sound speeds of the ionized and of the the neutral fluids are defined as c_{ie}=\sqrt{\gamma P_{ie}^{(0)}/\rho_{i}^{(0)}} and c_{S,n}=\sqrt{\gamma P_{n}^{(0)}/\rho_{n}^{(0)}}, respectively. In the fully ionized singlefluid case, the secondorder perturbations of pressure and density are related by the expression P_{ie}^{(2)}=c_{ie}^{2}\rho^{(2)} (see, e.g., Hollweg, 1971; Rankin et al., 1994). When multifluid plasmas are considered, that relation is not accurate because of the heat transfer terms in the evolution equation of pressure (see Equation (3) of Paper I). Nevertheless, for the purposes of this analytical study, it can be taken as a good approximation. Thus, assuming in the same way that P_{n}^{(2)}\approx c_{S,n}^{2}\rho_{n}^{(2)} and combining Equations (A9)(A12), it is possible to obtain the following equation that describes the secondorder perturbations of the density of ions (a similar equation can be cast for neutrals and for the xcomponent of the velocities of ions and neutrals):
\displaystyle\left[\frac{\partial^{4}}{\partial t^{4}}+\left(\nu_{in}+\nu_{ni}% \right)\frac{\partial^{3}}{\partial t^{3}}\left(c_{S,n}^{2}+c_{ie}^{2}\right)% \frac{\partial^{2}}{\partial t^{2}}\frac{\partial^{2}}{\partial x^{2}}\left(% \nu_{in}c_{S,n}^{2}+\nu_{ni}c_{ie}^{2}\right)\frac{\partial}{\partial t}\frac{% \partial^{2}}{\partial x^{2}}+c_{ie}^{2}c_{S,n}^{2}\frac{\partial^{4}}{% \partial x^{4}}\right]\rho_{i}^{(2)}=  
\displaystyle\left(\frac{\partial^{2}}{\partial t^{2}}\frac{\partial^{2}}{% \partial x^{2}}c_{S,n}^{2}\frac{\partial^{4}}{\partial x^{4}}+\nu_{ni}\frac{% \partial}{\partial t}\frac{\partial^{2}}{\partial x^{2}}\right)\left(\frac{B_{% \bot}^{2}}{2\mu_{0}}\right).  (A13) 
An interesting limiting case of the previous equation can be found if \nu_{ni} is assumed to tend to infinity, which corresponds to a strong coupling between the two fluids. The following expression is obtained:
\left[\left(1+\chi\right)\frac{\partial^{3}}{\partial t^{3}}\left(\chi c_{S,n% }^{2}+c_{ie}^{2}\right)\frac{\partial}{\partial t}\frac{\partial^{2}}{\partial x% ^{2}}\right]\rho_{i}^{(2)}=\frac{\partial}{\partial t}\frac{\partial^{2}}{% \partial x^{2}}\left(\frac{B_{\bot}^{2}}{2\mu_{0}}\right),  (A14) 
where the relation \nu_{in}/\nu_{ni}=\chi has been used. The integration with respect to time leads to the inhomogeneous wave equation, with a driving term in the righthand side, namely
\left(\frac{\partial^{2}}{\partial t^{2}}\widetilde{c}_{S}\frac{\partial^{2}}% {\partial x^{2}}\right)\rho_{i}^{(2)}=\frac{\partial^{2}}{\partial x^{2}}\left% (\frac{B_{\bot}^{2}}{2\mu_{0}(1+\chi)}\right),  (A15) 
where an integration constant has been taken equal to zero and \widetilde{c}_{S} is an effective sound speed given by
\widetilde{c}_{S}=\left(\frac{c_{ie}^{2}+\chi c_{S,n}^{2}}{1+\chi}\right)^{1/2}.  (A16) 
Following a similar procedure, the differential equation for the secondorder perturbation in the longitudinal component of the velocity is given by
\left(\frac{\partial^{2}}{\partial t^{2}}\widetilde{c}_{S}^{2}\frac{\partial^% {2}}{\partial x^{2}}\right)V_{i,x}^{(2)}=\frac{1}{\rho_{i}}\frac{\partial}{% \partial t}\frac{\partial}{\partial x}\left(\frac{B_{\bot}^{2}}{2\mu_{0}(1+% \chi)}\right),  (A17) 
From Equation (A13) it is also possible to recover the differential equation that describes the secondorder perturbations of density in a fully ionized plasma. If the collision frequencies are taken equal to zero (meaning that neutrals are decoupled and do not interact with ions), it is possible to rewrite Equation (A13) as
\left(\frac{\partial^{2}}{\partial x^{2}}c_{S,n}^{2}\frac{\partial^{2}}{% \partial x^{2}}\right)\left(\frac{\partial^{2}}{\partial x^{2}}c_{ie}^{2}% \frac{\partial^{2}}{\partial x^{2}}\right)\rho_{i}^{(2)}=\left(\frac{\partial^% {2}}{\partial x^{2}}c_{S,n}^{2}\frac{\partial^{2}}{\partial x^{2}}\right)% \frac{\partial^{2}}{\partial x^{2}}\left(\frac{B_{\bot}^{2}}{2\mu_{0}}\right),  (A18) 
which leads to
\left(\frac{\partial^{2}}{\partial t^{2}}c_{ie}^{2}\frac{\partial^{2}}{% \partial x^{2}}\right)\rho_{i}^{(2)}=\frac{\partial^{2}}{\partial x^{2}}\left(% \frac{B_{\bot}^{2}}{2\mu_{0}}\right),  (A19) 
an equation that has already been derived by Hollweg (1971), Tikhonchuk et al. (1995) or Terradas & Ofman (2004). It can be seen that Equations (A15) and (A19) represent the same type of behavior, with differences appearing in the velocity of propagation of waves and the amplitude of the driving term. These are two effects caused by the ionneutral interaction.
A.1 Standing waves
In this section, the properties of nonlinear standing waves in a twofluid partially ionized plasma are analyzed.
If the initial perturbation applied to the equilibrium state is given by
V_{y}^{(1)}(x,t=0)=V_{y,0}\cos(k_{x}x),  (A20) 
and the strongly coupled limit is applied to Equations (A7) and (A8), the firstorder perturbation of the magnetic field is
B_{\bot}(x,t)=\frac{B_{0}}{\widetilde{c}_{\mathrm{A}}}V_{y,0}\sin(\widetilde{% c}_{\mathrm{A}}k_{x}t)\sin(k_{x}x),  (A21) 
with \widetilde{c}_{\mathrm{A}} the Alfvén speed modified by the inclusion of the inertia of neutrals, i.e., \widetilde{c}_{\mathrm{A}}=B_{0}/\sqrt{\mu_{0}\rho_{i,0}(1+\chi)}. Then, using the initial conditions \rho_{i}^{(2)}(x,t=0)=0 and \frac{\partial}{\partial t}\rho_{i}^{(2)}(x,t=0)=0, respectively, the solution to Equation (A15) can be computed as
\rho_{i}^{(2)}(x,t)=\frac{1}{2\widetilde{c}_{S}}\int_{0}^{t}\int_{x\widetilde% {c}_{S}(t\tau)}^{x+\widetilde{c}_{S}(t\tau)}\frac{\partial^{2}}{\partial x^{% 2}}\left[\frac{B_{\bot}^{2}(\xi,\tau)}{2\mu_{0}(1+\chi)}\right]d\xi d\tau.  (A22) 
The only speed that explicitly appears in Equation (A22) is the effective sound speed, \widetilde{c}_{S}. However, since the driving wave is assumed to be Alfvénic, B_{\bot}^{2} is a function of the Alfvén speed. Hence, \rho_{i}^{(2)}(x,t) also depends on that speed. Finally, the secondorder perturbation of the density of ions is given by
\rho_{i}^{(2)}(x,t)=\frac{B_{0}^{2}V_{y,0}^{2}\left[\widetilde{c}_{\mathrm{A}}% ^{2}\bar{c}_{S}^{2}+\bar{c}_{S}^{2}\cos(2\widetilde{c}_{\mathrm{A}}k_{x}t)% \widetilde{c}_{\mathrm{A}}^{2}\cos(2\bar{c}_{S}k_{x}t)\right]\cos(2k_{x}x)}{8% \widetilde{c}_{\mathrm{A}}^{2}\bar{c}_{S}^{2}(\widetilde{c}_{\mathrm{A}}^{2}% \bar{c}_{S}^{2})\mu_{0}(1+\chi)}.  (A23) 
The resulting perturbation is the combination of two standing modes with frequencies 2\widetilde{c}_{\mathrm{A}}k_{x} and 2\widetilde{c}_{S}k_{x}, respectively, and whose wavenumber is twice the wavenumber of the original perturbation. The solution for the fully ionized case is recovered by substituting \widetilde{c}_{S} with c_{ie}, \widetilde{c}_{\mathrm{A}} with c_{\mathrm{A}} and taking \chi=0.
If the sound speed is much lower than the Alfvén speed, as it occurs in the simulations performed in this work, Equation (A23) can be approximated as
\rho_{i}^{(2)}(x,t)\approx\frac{B_{0}^{2}V_{y,0}^{2}}{8\mu_{0}\widetilde{c}_{% \mathrm{A}}^{2}(1+\chi)}\left[\frac{1\cos(2\widetilde{c}_{S}k_{x}t)}{% \widetilde{c}_{S}^{2}}\right]\cos(2k_{x}x),  (A24) 
which shows that the perturbation is dominated by the oscillation mode associated with the weighted sound speed.
Then, the relative variation of density, which in Figures 13 is represented as \Delta\rho/\rho_{0}, can be computed as the ratio between the secondorder perturbation and the background density. Hence,
\frac{\Delta\rho_{i}}{\rho_{i,0}}\equiv\frac{\rho_{i}^{(2)}(x,t)}{\rho_{i,0}}% \approx\frac{V_{y,0}^{2}}{8\widetilde{c}_{S}^{2}}\left[1\cos(2\widetilde{c}_{% S}k_{x}t)\right]\cos(2k_{x}x),  (A25) 
An interesting conclusion can be extracted from the previous equation: since the relative variation of density is proportional to V_{y,0}^{2}/\widetilde{c}_{S}^{2} for partially ionized plasmas while it is proportional to V_{y,0}^{2}/c_{ie}^{2} for fully ionized fluids and c_{ie}>\widetilde{c}_{S}, the relative variation of density is larger when the effect of partial ionization is taken into account. This is an important result caused by partial ionization.
If the wavenumber of the perturbation increases, the frequency of the Alfvén wave increases as well and the coupling between the two fluids is not as strong as for smaller wavenumbers. Hence, it would be expected that Equation (A25) becomes inaccurate at larger wavenumbers. Moreover, it has been shown in Papers I and II that Hall’s term should be taken into account in the large wavenumber range.
The threefluid simulations represented in Figures 2 and 3 show that, under the chosen physical parameters, the friction due to ionneutral collisions is more efficient in attenuating the Alfvénic waves than the acoustic modes. For instance, it can be checked that in Figure 3 the firstorder Alfvén wave has almost disappeared after t=0.5\ \mathrm{s}, but the secondorder perturbation in the xcomponent of the velocity lasts for a longer time. In a twofluid plasma, the oscillation frequency and damping rate of the remaining secondorder wave may be obtained from Equation (A13) in the following way. Since the driving wave, i.e., the firstorder Alfvén wave, vanishes due to collisions, after a given time the term on the righthand side of Equation (A13) becomes equal to zero. Then, the remaining oscillations are governed by the homogeneous version of the differential equation, with the initial conditions given by the wave previously induced by the driver. After the primary Alfvén wave is completely damped, the secondorder perturbation of the density of ions can be expressed as
\rho_{i}^{(2)}\sim\exp\left[i(\omega t+\kappa x)\right],  (A26) 
where, in this case, the wavenumber is twice the wavenumber of the original driving wave, i.e., \kappa=2k_{x}. This procedure leads to the following dispersion relation,
\omega^{4}+i\left(\nu_{ni}+\nu_{in}\right)\omega^{3}\kappa^{2}\left(c_{S,n}^{% 2}+c_{ie}^{2}\right)\omega^{2}i\kappa^{2}\left(\nu_{in}c_{S,n}^{2}+\nu_{ni}c_% {ie}^{2}\right)\omega+c_{ie}^{2}c_{S,n}^{2}\kappa^{4}=0,  (A27) 
which depends on the sound speeds but not on the Alfvén speed. This is the same dispersion relation that would be obtained for linear acoustic waves in a twospecies fluid in which only the collisional interaction between ions and neutrals is taken into account and the influence of magnetic fields is neglected (see, e.g., Vranjes & Poedts, 2010). It coincides with Equation (9) from Vranjes & Poedts (2010) if the factors proportional to the electronneutral collision frequency of that formula are neglected, and it can also be recovered from Equation (47) of Soler et al. (2013a), where magnetoacoustic waves in partially ionized plasmas have been studied, if the Alfvén speed is set equal to zero.
It must be noted that for a certain range of collision frequencies, the driving wave may last more than the acoustic wave and, strictly, the dispersion relation, Equation (A27) should not be applicable because the driver is still working. This is a consequence of the damping due to ionneutral collisions being most efficient when the oscillation frequency is similar to the collision frequency (Zaqarashvili et al., 2011; Soler et al., 2013b). Since \widetilde{c}_{S}\ll\widetilde{c}_{\mathrm{A}}, the acoustic modes are more damped than the Alfvénic ones at low collision frequencies and the opposite would occur at high frequencies. Nevertheless, as shown by Equations (A23) and (A24), if the Alfvén speed is much larger than the sound speed, the secondorder oscillation is dominated by the acoustic mode. Hence, the results from Equation (A27) are still good approximations at any range of collision frequencies.
A.2 Propagating waves
Here, we study the case of nonlinear propagating waves. A linearly polarized Alfvén wave is driven at the boundary x=0 by applying the following boundary condition to the ycomponent of the velocity:
V_{y}^{(1)}(0,t)=V_{y,0}\sin\left(\omega t\right),  (A28) 
where V_{y,0} is proportional to the Alfvén speed, c_{\mathrm{A}}.
If we consider the case of strong coupling between ions and neutrals, i.e., \nu_{ni}\to\infty, and focus only on waves that propagate along the positive direction of the xaxis, Equations (A7) and (A8) give the following firstorder perturbations of the velocity and magnetic field:
V_{y}^{(1)}(x,t)=\left\{\begin{array}[]{l c}V_{y,0}\sin\left[\omega\left(t% \frac{x}{\widetilde{c}_{\mathrm{A}}}\right)\right],&\text{if}\ x\leq\widetilde% {c}_{\mathrm{A}}t,\\ 0,&\text{if}\ x>\widetilde{c}_{\mathrm{A}}t,\end{array}\right.  (A29) 
B_{\bot}(x,t)=\left\{\begin{array}[]{l c}\frac{B_{0}V_{y,0}}{\widetilde{c}_{% \mathrm{A}}}\sin\left[\omega\left(t\frac{x}{\widetilde{c}_{\mathrm{A}}}\right% )\right],&\text{if $x\leq\widetilde{c}_{\mathrm{A}}t$},\\ 0,&\text{if $x>\widetilde{c}_{\mathrm{A}}t$}.\end{array}\right.  (A30) 
The expression for the magnetic field perturbation can be inserted into Equation (A17) to obtain the secondorder perturbation for the longitudinal velocity. This perturbation is the combination of the solutions of the homogeneous and inhomogeneous wave equation. The inhomogeneous solution can be obtained by assuming that it is given by V_{I}^{(2)}(x,t)=A_{1}\cos\left(\delta_{1}\omega t\delta_{2}x\right)+A_{2} and inserting this expression into Equation (A17). This step allows to obtain the values of A_{1}, \delta_{1}, and \delta_{2}. The value of A_{2} is computed by imposing that V_{NH}^{(2)}(0,0)=0. Then, the inhomogeneous solution is given by
V_{I}^{(2)}(x,t)=\frac{V_{y,0}^{2}\widetilde{c}_{\mathrm{A}}}{4\left(% \widetilde{c}_{\mathrm{A}}^{2}\widetilde{c}_{S}^{2}\right)}\left\{1\cos\left% [2\omega\left(t\frac{x}{\widetilde{c}_{\mathrm{A}}}\right)\right]\right\}.  (A31) 
The homogeneous solution describes a wave that propagates at the effective sound speed, \widetilde{c}_{S}. It is obtained by assuming that V_{H}^{(2)}(x,t)=A_{3}\cos\left[\delta_{3}\omega\left(tx/\widetilde{c}_{S}% \right)\right]+A_{4} and imposing that the complete solution fulfills the boundary condition V_{i,x}^{(2)}(0,t)=0, which means that there is no mass inflow. Hence,
V_{H}^{(2)}(x,t)=\frac{V_{y,0}^{2}\widetilde{c}_{\mathrm{A}}}{4\left(% \widetilde{c}_{\mathrm{A}}^{2}\widetilde{c}_{S}^{2}\right)}\left\{\cos\left[2% \omega\left(t\frac{x}{\widetilde{c}_{S}}\right)\right]1\right\}.  (A32) 
Under the physical conditions analyzed in this work, Alfvén waves propagate faster than sound waves, as shown by Figure 11. Consequently, the secondorder perturbation of the longitudinal component of the velocity is given by
V_{i,x}^{(2)}(x,t)=\left\{\begin{array}[]{l c}V_{H}^{(2)}(x,t)+V_{I}^{(2)}(x,t% ),&\text{if $x\leq\widetilde{c}_{S}t$},\\ V_{I}^{(2)}(x,t),&\text{if $\widetilde{c}_{S}t<x<\widetilde{c}_{\mathrm{A}}t$}% ,\\ 0,&\text{if $x>\widetilde{c}_{\mathrm{A}}t$}.\end{array}\right.  (A33) 
The secondorder perturbation of density can be computed through the continuity equation, i.e., Equation (A10), taking into account that \rho_{i}^{(2)}(0,0)=0. This leads to
\rho_{i}^{(2)}(x,t)=\left\{\begin{array}[]{l c}\rho_{H}^{(2)}(x,t)+\rho_{I}^{(% 2)}(x,t),&\text{if $x\leq\widetilde{c}_{S}t$},\\ \rho_{I}^{(2)}(x,t),&\text{if $\widetilde{c}_{S}t<x<\widetilde{c}_{\mathrm{A}}% t$},\\ 0,&\text{if $x>\widetilde{c}_{\mathrm{A}}t$},\end{array}\right.  (A34) 
where
\rho_{I}^{(2)}(x,t)=\frac{\rho_{i,0}V_{y,0}^{2}}{4\left(\widetilde{c}_{\mathrm% {A}}^{2}\widetilde{c}_{S}^{2}\right)}\left\{1\cos\left[2\omega\left(t\frac{% x}{\widetilde{c}_{\mathrm{A}}}\right)\right]\right\}  (A35) 
and
\rho_{H}^{(2)}(x,t)=\frac{\rho_{i,0}V_{y,0}^{2}\widetilde{c}_{\mathrm{A}}}{4% \widetilde{c}_{S}\left(\widetilde{c}_{\mathrm{A}}^{2}\widetilde{c}_{S}^{2}% \right)}\left\{\cos\left[2\omega\left(t\frac{x}{\widetilde{c}_{S}}\right)% \right]1\right\}.  (A36) 
The solutions for the fully ionized case are recovered by substituting \widetilde{c}_{\mathrm{A}} and \widetilde{c}_{S} with c_{\mathrm{A}} and c_{ie}, respectively. In that case, Equations (A31) and (A32) coincide with Equations (25) and (26) from Zheng et al. (2016), who studied the propagation of nonlinear Alfvén waves in ideal MHD plasmas.
Equations (A31), (A32), (A35), and (A36) show that the oscillation frequency of the secondorder perturbations is twice the frequency of the driver and the amplitude has a quadratic dependence on the amplitude of the driving wave. Since \widetilde{c}_{S}<\widetilde{c}_{\mathrm{A}}, the amplitude of \rho_{H}^{(2)} is larger than that of \rho_{I}^{(2)}. In addition, it can be checked that \rho_{H}^{(2)} is lower or equal to zero, while \rho_{I}^{(2)} is larger or equal to zero. The net effect of the combination of these two propagating perturbations is that the total amount of matter in the region x<\widetilde{c}_{S}t decreases.
References
 Arber et al. (2016) Arber, T. D., Brady, C. S., & Shelyag, S. 2016, ApJ, 817, 94
 Barnes & Hollweg (1974) Barnes, A., & Hollweg, J. V. 1974, JGR, 79, 2302
 Belcher & Davis (1971) Belcher, J. W., & Davis, Jr., L. 1971, JGR, 76, 3534
 Belcher et al. (1969) Belcher, J. W., Davis, Jr., L., & Smith, E. J. 1969, JGR, 74, 2302
 Bona et al. (2009) Bona, C., BonaCasas, C., & Terradas, J. 2009, Journal of Computational Physics, 228, 2266
 Cohen & Kulsrud (1974) Cohen, R. H., & Kulsrud, R. M. 1974, Physics of Fluids, 17, 2215
 Coleman (1968) Coleman, Jr., P. J. 1968, ApJ, 153, 371
 Esser et al. (1986) Esser, R., Habbal, S. R., Withbroe, G. L., & Leer, E. 1986, JGR, 91, 2950
 Eto et al. (2002) Eto, S., Isobe, H., Narukage, N., et al. 2002, PASJ, 54, 481
 Gilbert (2015) Gilbert, H. 2015, in Astrophysics and Space Science Library, Vol. 415, Solar Prominences, ed. J.C. Vial & O. Engvold, 157
 Gilbert et al. (2008) Gilbert, H. R., Daou, A. G., Young, D., Tripathi, D., & Alexander, D. 2008, ApJ, 685, 629
 Goodman (2011) Goodman, M. L. 2011, ApJ, 735, 45
 Haerendel (1992) Haerendel, G. 1992, Nature, 360, 241
 Heinzel (2015) Heinzel, P. 2015, in Astrophysics and Space Science Library, Vol. 415, Solar Prominences, ed. J.C. Vial & O. Engvold, 103
 Heinzel et al. (2010) Heinzel, P., Anzer, U., & Gunár, S. 2010, Mem. Soc. Astron. Italiana, 81, 654
 Heinzel et al. (2015) Heinzel, P., Gunár, S., & Anzer, U. 2015, A&A, 579, A16
 Hollweg (1971) Hollweg, J. V. 1971, JGR, 76, 5155
 Jing et al. (2003) Jing, J., Lee, J., Spirock, T. J., et al. 2003, Astrophysical Journal Letters, 584, L103
 Kamaya & Nishi (1998) Kamaya, H., & Nishi, R. 1998, ApJ, 500, 257
 Karpen et al. (2017) Karpen, J. T., DeVore, C. R., Antiochos, S. K., & Pariat, E. 2017, ApJ, 834, 62
 Khomenko & Collados (2012) Khomenko, E., & Collados, M. 2012, ApJ, 747, 87
 Kulsrud & Pearce (1969) Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445
 Kumar & Roberts (2003) Kumar, N., & Roberts, B. 2003, Sol. Phys., 214, 241
 Leake et al. (2005) Leake, J. E., Arber, T. D., & Khodachenko, M. L. 2005, A&A, 442, 1091
 Lou (1993) Lou, Y.Q. 1993, JGR, 98, 3563
 Luna & Karpen (2012) Luna, M., & Karpen, J. 2012, Astrophysical Journal Letters, 750, L1
 Martin et al. (1997) Martin, C. E., Heyvaerts, J., & Priest, E. R. 1997, A&A, 326, 1176
 MartínezGómez et al. (2016) MartínezGómez, D., Soler, R., & Terradas, J. 2016, ApJ, 832, 101,
 MartínezGómez et al. (2017) —. 2017, ApJ, 837, 80,
 Matsumoto & Shibata (2010) Matsumoto, T., & Shibata, K. 2010, ApJ, 710, 1857
 McIntosh et al. (2011) McIntosh, S. W., de Pontieu, B., Carlsson, M., et al. 2011, Nature, 475, 477
 Moreton & Ramsey (1960) Moreton, G. E., & Ramsey, H. E. 1960, PASP, 72, 357
 Moses et al. (1997) Moses, D., Clette, F., Delaboudinière, J.P., et al. 1997, SoPh, 175, 571
 Mouschovias et al. (2011) Mouschovias, T. C., Ciolek, G. E., & Morton, S. A. 2011, MNRAS, 415, 1751
 Murawski & Roberts (1993) Murawski, K., & Roberts, B. 1993, SoPh, 145, 65
 Ofman & Davila (1998) Ofman, L., & Davila, J. M. 1998, JGR, 103, 23677
 Okamoto et al. (2004) Okamoto, T. J., Nakai, H., Keiyama, A., et al. 2004, ApJ, 608, 1124
 Oliver et al. (1998) Oliver, R., Murawski, K., & Ballester, J. L. 1998, A&A, 330, 726
 Pandey & Wardle (2008) Pandey, B. P., & Wardle, M. 2008, MNRAS, 385, 2269
 Parenti et al. (2006) Parenti, S., Buchlin, E., Cargill, P. J., Galtier, S., & Vial, J.C. 2006, ApJ, 651, 1219
 Parenti & Vial (2007) Parenti, S., & Vial, J.C. 2007, A&A, 469, 1109
 Parnell & De Moortel (2012) Parnell, C. E., & De Moortel, I. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 3217
 Piddington (1956) Piddington, J. H. 1956, MNRAS, 116, 314
 Pudritz (1990) Pudritz, R. E. 1990, ApJ, 350, 195
 Ramsey & Smith (1966) Ramsey, H. E., & Smith, S. F. 1966, AJ, 71, 197
 Rankin et al. (1994) Rankin, R., Frycz, P., Tikhonchuk, V. T., & Samson, J. C. 1994, JGR, 99, 21291
 Ruderman (2006) Ruderman, M. S. 2006, Philosophical Transactions of the Royal Society of London Series A, 364, 485
 Schunk (1977) Schunk, R. W. 1977, Reviews of Geophysics and Space Physics, 15, 429
 Soler et al. (2013a) Soler, R., Carbonell, M., & Ballester, J. L. 2013a, ApJs, 209, 16
 Soler et al. (2013b) Soler, R., Carbonell, M., Ballester, J. L., & Terradas, J. 2013b, ApJ, 767, 171
 Soler et al. (2016) Soler, R., Terradas, J., Oliver, R., & Ballester, J. L. 2016, A&A, 592, A28
 Soler et al. (2017) —. 2017, ApJ, 840, 20
 Song & Vasyliūnas (2011) Song, P., & Vasyliūnas, V. M. 2011, Journal of Geophysical Research (Space Physics), 116, A09104
 Spitzer (1956) Spitzer, L. 1956, Physics of Fully Ionized Gases
 Srivastava et al. (2017) Srivastava, A. K., Shetye, J., Murawski, K., et al. 2017, Scientific Reports, 7, 43147
 Suzuki (2008) Suzuki, T. K. 2008, Nonlinear Processes in Geophysics, 15, 295
 Suzuki (2011) —. 2011, SSRv, 158, 339
 Terradas & Ofman (2004) Terradas, J., & Ofman, L. 2004, ApJ, 610, 523
 Thompson et al. (1998) Thompson, B. J., Plunkett, S. P., Gurman, J. B., et al. 1998, Geophys. Res. Lett., 25, 2465
 Tikhonchuk et al. (1995) Tikhonchuk, V. T., Rankin, R., Frycz, P., & Samson, J. C. 1995, Physics of Plasmas, 2, 501
 Tomczyk et al. (2007) Tomczyk, S., McIntosh, S. W., Keil, S. L., et al. 2007, Science, 317, 1192
 Tu & Song (2013) Tu, J., & Song, P. 2013, ApJ, 777, 53
 Verwichte et al. (1999) Verwichte, E., Nakariakov, V. M., & Longbottom, A. W. 1999, Journal of Plasma Physics, 62, 219
 Voitenko & Goossens (2002) Voitenko, Y., & Goossens, M. 2002, SoPh, 209, 37
 Vranjes & Poedts (2010) Vranjes, J., & Poedts, S. 2010, Physics of Plasmas, 17, 022104
 Watanabe (1961) Watanabe, T. 1961, Canadian Journal of Physics, 39, 1197
 Zaqarashvili et al. (2011) Zaqarashvili, T. V., Khodachenko, M. L., & Rucker, H. O. 2011, A&A, 529, A82
 Zheng et al. (2016) Zheng, J., Chen, Y., & Yu, M. 2016, Phys. Scr, 91, 015601