Multi-fluid Approach to High-frequency Waves in Plasmas. III. Nonlinear Regime and Plasma Heating

Alfvén waves are usually classified into two categories depending on their velocity amplitudes. Small-amplitude waves have velocity amplitudes that are much smaller than the Alfvén speed. Conversely, the velocity amplitudes of large-amplitude waves are not negligible in comparison with the Alfvén speed. In Martínez-Gómez et al. (2016, 2017; hereafter, Papers I and II, respectively), we studied the properties of small-amplitude waves in fully and partially ionized plasmas of the solar atmosphere by considering the linear regime of a multi-fluid 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 large-amplitude perturbations are considered.

Solar prominences and the solar wind are examples of solar plasmas in which large-amplitude waves have been detected. Large-amplitude 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 over a few periods. These large-amplitude oscillations are rare events, but in the last few years a growing number of observations have 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) have detected large-amplitude 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), and 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), and 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 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; and Karpen et al. (2017), whose results suggest that coronal-hole jets are possible origins of 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 higher-order perturbations. The series is truncated at some given order, and systems of equations are derived for the perturbations, while higher-order effects are ignored. This procedure was followed, e.g., by Hollweg (1971), who studied second-order effects of Alfvén waves, and by Rankin et al. (1994), Tikhonchuk et al. (1995), and Verwichte et al. (1999), who examined the properties of up to third-order perturbations. These 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, third-order 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 single fluid. The assumption of full 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 single-fluid models is only appropriate when the phenomenon of interest is associated with low frequencies, i.e., much lower than the ion cyclotron frequencies in fully ionized plasmas or the ion–neutral collision frequencies in partially ionized plasmas. Conversely, at higher frequencies, multi-fluid approaches are required due to the fact that the components of the plasma are not strongly coupled and react to perturbations on different timescales.

In the present work, the multi-fluid model described in Paper I is applied to the investigation of nonlinear waves in partially ionized plasmas, and special attention is paid to the heating due to ion–neutral collisions. The issue of heating is of great interest in the solar atmosphere (see, e.g., Goodman 2011; Song & Vasyliūnas 2011; Khomenko & Collados 2012; Parnell & De Moortel 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 of 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 ion–neutral collisional interaction is one of the possible mechanisms. The topic of heating by means of ion–neutral collisions was briefly examined in Paper II when small-amplitude 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 large-amplitude waves are considered. The model used here also considers Coulomb collisions, magnetic diffusivity, and the effects of Hall's current, and accounts for the cyclotron motion of ions and has been shown to be greatly relevant 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, large-amplitude impulsive perturbations are considered, and the heating due to ion–neutral 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 two-fluid plasmas.

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 five-moment 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 the magnetic resistivity or Ohm's diffusion. Interested readers are referred to Paper I for the detailed derivation and discussion of the governing multi-fluid 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 km over the photosphere and with a gas pressure of P_g = 0.005 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, H, and He, respectively. Ionized helium has a residual abundance at the cool temperatures of the prominence cores. So, for simplicity, we assume helium to be fully neutral. The presence of ionized helium is, however, important for prominence-to-corona 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}\,{{\rm{m}}}^{-3}$, ${n}_{{\rm{H}}}=2\times {10}^{16}\,{{\rm{m}}}^{-3}$, ${n}_{\mathrm{He}}=2\times {10}^{15}\,{{\rm{m}}}^{-3}$, and T₀ = 10,000 K, respectively. The collision frequencies of protons with neutral hydrogen, of protons with neutral helium, and of hydrogen with helium are ${\nu }_{p{\rm{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 }_{\mathrm{st}}$, between two species s and t, and the friction coefficient, α_st, are related by ${\nu }_{\mathrm{st}}={\alpha }_{\mathrm{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, ${{\boldsymbol{B}}}_{0}$, along the x-direction is considered. A typical value of the magnetic field strength in quiescent prominences is B₀ = 10 G. The fundamental standing mode of the transverse Alfvén waves is excited by applying the initial perturbation

$$\begin{eqnarray}&&{V}_{s,y}(x,t=0)={V}_{y,0}\cos ({k}_{x}x)\end{eqnarray} \tag{ 1 }$$

to every species s of the plasma, where ${V}_{s,y}$ is the y-component 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,{\boldsymbol{B}}(x,t=0)=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}\,{\rm{m}}$. Thus, the wavenumber of the fundamental mode is ${k}_{x}=\pi /(2l)=\pi /(5\times {10}^{5})\,{{\rm{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 the variables are extrapolated. The amplitude of the perturbation is given by ${V}_{y,0}=2.5\times {10}^{-2}{c}_{{\rm{A}}}$, where ${c}_{{\rm{A}}}=| {{\boldsymbol{B}}}_{0}| /\sqrt{{\mu }_{0}{\rho }_{p}}$ is the Alfvén speed, with μ₀ 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}_{{\rm{A}}}\approx 184\,\mathrm{km}\,{{\rm{s}}}^{-1}$.

The top row of Figure 1 shows several snapshots of the evolution of the y-component 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 y-component of the magnetic field and on the z-components 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 entire 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 these reasons, to analyze the properties of the transverse waves, we focus only on the y-component 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 illustrated in the top-left panel of Figure 2. This is the same behavior as the one already explained in Paper II for small-amplitude 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 x-component of the velocity, normalized with respect to the amplitude of the driver, ${V}_{y,0}$, at various time steps. Aside 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 (2{k}_{x}x)$. Furthermore, the top-right panel of Figure 2 shows that the oscillation in V_x does not attenuate as fast as the oscillation in V_y.

Zoom In Zoom Out Reset image size

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., ${\rm{\Delta }}\rho /{\rho }_{0}$, with ${\rm{\Delta }}\rho \equiv \rho (x,t)-{\rho }_{0}$, and the ratio between the perturbation of the temperature and the initial temperature, ${\rm{\Delta }}T/{T}_{0}$, with ${\rm{\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}_{{\rm{B}}}{T}_{s}$, where ${k}_{{\rm{B}}}$ is Boltzmann's constant. It can be checked that both ${\rm{\Delta }}\rho /{\rho }_{0}$ and ${\rm{\Delta }}T/{T}_{0}$ are proportional to $\cos (2{k}_{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 this process is later reversed and an oscillation appears. The amplitude of this variation of the 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 ion–neutral 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 the 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 wave of frequency ω, we find that $\omega \approx 0.67\,\mathrm{rad}\,{{\rm{s}}}^{-1}$ (which corresponds to a period of 9.4 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π (it is common practice to compare ω directly with ν_st but, rigorously speaking, this comparison is not correct because those quantities are expressed in different units), we find that ω/(2π) < ν_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 x-component 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}\,{{\rm{s}}}^{-1}$ and ${\omega }_{2}\approx 1.34\,\mathrm{rad}\,{{\rm{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 ${\rm{\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 and reaches a maximum; the fluctuation can then 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 ${\rm{\Delta }}\rho /{\rho }_{0}$ and V_x. Finally, panel (d) shows a rising 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 ion–neutral 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 )\gt {\nu }_{p\mathrm{He}}$. To that end, we perform a simulation with a larger wavenumber, ${k}_{x}=\pi /(5\times {10}^{3})\,{{\rm{m}}}^{-1}$. The dispersion relation (see Paper II) predicts a frequency $\omega \approx 75.14\,\mathrm{rad}\,{{\rm{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{\rm{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 x-component of the 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 quickly. 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 multi-fluid effect that cannot be captured with the usual single-fluid treatments.

Zoom In Zoom Out Reset image size

Figure 3(c) shows that the density only increases at the center of the domain for a very short time. Then, the relative variation of density becomes negative and oscillates about Δρ/ρ₀ ≈ −0.02. Hence, the net result of this nonlinear effect is that matter is displaced from the central part of the domain and directed toward 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 the temperature and pressure is caused by the fast dissipation of the Alfvén wave due to the collisional friction. This issue will be addressed in more detail later. It is also remarkable that during the first steps of the simulation, neutral helium reaches a higher 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 on 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 multi-fluid equations. However, a three-fluid 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 a simpler scenario can be found in the Appendix. Here, we just mention its main conclusions:

• 1.  
  A standing Alfvén wave nonlinearly generates two second-order perturbations in the density and pressure, as well as in the longitudinal component of the velocity. The wavenumber of those perturbations, κ, is twice the wavenumber of the original wave, i.e., $\kappa =2{k}_{x};$
• 2.  
  The frequencies of the second-order perturbations are given by $2{\widetilde{c}}_{S}{k}_{x}$ and $2{\widetilde{c}}_{{\rm{A}}}{k}_{x}$, where ${\widetilde{c}}_{S}$ and ${\widetilde{c}}_{{\rm{A}}}$ are the effective sound speed and the modified Alfvén speed, respectively. The former is given by the positive square root of

  $$\begin{eqnarray}&&{\widetilde{c}}_{S}^{2}=\displaystyle \frac{{\displaystyle \sum }_{t}{\rho }_{t}{c}_{S,t}^{2}}{{\displaystyle \sum }_{t}{\rho }_{t}},\end{eqnarray} \tag{ 2 }$$

  where ${c}_{S,t}=\sqrt{\gamma {P}_{t,0}/{\rho }_{t,0}}$ is the sound speed of species t, with γ the adiabatic constant and ${P}_{t,0}$ the equilibrium value of pressure. The modified Alfvén speed is given by

  $$\begin{eqnarray}&&{\widetilde{{\boldsymbol{c}}}}_{{\boldsymbol{A}}}=\displaystyle \frac{{{\boldsymbol{c}}}_{{\rm{A}}}}{\sqrt{1+{\displaystyle \sum }_{{\boldsymbol{t}}\ne {\boldsymbol{p}}}{{\boldsymbol{\chi }}}_{{\boldsymbol{t}}}}},\end{eqnarray} \tag{ 3 }$$

  where ${\chi }_{t}={\rho }_{t}/{\rho }_{p};$
• 3.  
  And, if ${\widetilde{c}}_{S}^{2}\ll {\widetilde{c}}_{{\rm{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 1–3. However, it must be noted that they correspond to a case in which the coupling between all of 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 second-order perturbations on the wavenumber. The results of this study are represented in Figure 4, where the normalized oscillation frequency of waves, ${\omega }_{R}/{\omega }_{{\rm{A}}}$, is plotted as a function of the wavenumber k_x. The normalization constant is the Alfvén frequency, ${\omega }_{{\rm{A}}}={k}_{x}{c}_{{\rm{A}}}$.

Zoom In Zoom Out Reset image size

The black lines in 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 R-mode and the dashed line represents the L-mode. 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 first-order perturbations.

Regarding the second-order modes, their frequency, denoted by ${\omega }_{\mathrm{sim}}^{(2)}$, has three clearly different regimes depending on the wavenumber of the first-order perturbation. At small wavenumbers, the frequency is approximately given by ${\omega }_{\mathrm{sim}}^{(2)}\approx 2{\widetilde{c}}_{S}{k}_{x}$ (see Equation (40)). However, as k_x is increased, ω⁽²⁾ departs from that value. To understand this behavior, we take into account that at small wavenumbers, the oscillation frequency of the first-order perturbation, ${\omega }_{\mathrm{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 }_{\mathrm{sim}}^{(2)}/{\omega }_{{\rm{A}}}=2{\widetilde{c}}_{S}{k}_{x}/{\omega }_{{\rm{A}}}\approx 0.14$. At intermediate wavenumbers, ${\omega }_{\mathrm{sim}}^{(1)}$ is larger than ${\nu }_{p\mathrm{He}}$ and ${\nu }_{\mathrm{HHe}}$, but smaller than ${\nu }_{p{\rm{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 }_{\mathrm{sim}}^{(2)}/{\omega }_{{\rm{A}}}\approx 0.15$. Finally, at large wavenumbers, ${\omega }_{\mathrm{sim}}^{(1)}\gg {\nu }_{p{\rm{H}}}$ and the coupling between protons and hydrogen is weak. The sound speed of the proton fluid is c_ie and ${\omega }_{\mathrm{sim}}^{(2)}/{\omega }_{{\rm{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 $2{k}_{x}{c}_{S,{\rm{H}}}/{\omega }_{{\rm{A}}}\,\approx 0.13$ and $2{k}_{x}{c}_{S,\mathrm{He}}\approx 0.065$, respectively.

Furthermore, the simulations represented in Figures 1–3 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, ${\rm{\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 second-order 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 the magnetic pressure moves the plasma toward the nodes of the first-order magnetic field perturbation. The gradient of the second-order 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 timescales 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 the magnetic pressure has a shorter timescale 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 ion–neutral 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 second-order effects are relevant for the dynamics. However, if the amplitudes are increased, higher-order terms may also be of great importance. As detailed by Tikhonchuk et al. (1995), the higher-order 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 higher-order 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 x-component of the velocity at the point x = −l/2 and the variation of density at x = 0, respectively. The top left panel represents the first-order 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.15{c}_{{\rm{A}}}$ than for ${V}_{y,0}=0.025{c}_{{\rm{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.

Zoom In Zoom Out Reset image size

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 spatially average temperature, given by $1/(2l){\int }_{-l}^{l}T(x)\,{dx}$, is shown 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 second-order acoustic waves to heating is negligible. This is due to ion–neutral 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 }_{\mathrm{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.025{c}_{{\rm{A}}}$, the temperature rises up to ∼10,360 K (i.e., the variation is ΔT ≈ 360 K). For the amplitudes ${V}_{y,0}=0.05{c}_{{\rm{A}}}$, ${V}_{y,0}=0.1{c}_{{\rm{A}}}$, and ${V}_{y,0}=0.15{c}_{{\rm{A}}}$, the final temperatures are $11,470\,{\rm{K}}$ (ΔT ≈ 1470 K), 16,170 K (ΔT ≈ 6170 K), and 24,470 K (ΔT ≈ 14,470 K), respectively. Hence, the dependence of the increment of temperature on the amplitude of the perturbation is approximately quadratic.

Zoom In Zoom Out Reset image size

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 higher-order 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 second-order 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

$$\begin{eqnarray}&&{f}^{(1)}(x,t=0)\sim \exp \left[-{\left(\displaystyle \frac{x-{x}_{0}}{\sqrt{2}{\sigma }_{x}}\right)}^{2}\right],\end{eqnarray} \tag{ 4 }$$

where σ_x is the root-mean-square width and is related to the FWHM of the Gaussian through the formula $\mathrm{FWHM}\,=2\sqrt{2\mathrm{ln}2}{\sigma }_{x}$, and x₀ 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 y-component of the velocity of all species. The amplitude of the perturbation is ${V}_{y,0}=5\times {10}^{-2}{c}_{{\rm{A}}},$ and its width is $\mathrm{FWHM}=2\times {10}^{5}\,{\rm{m}}$. As expected, the initial pulse splits into two smaller Alfvénic pulses, with half the height of the initial pulse, and they propagate into opposite directions. There is 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}}_{{\rm{A}}}$. Notwithstanding this, 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.

Zoom In Zoom Out Reset image size

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 x-component 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 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}}_{{\rm{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 the 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 is later displaced from that location.

The third and fourth rows of Figure 8 represent the second-order perturbations of pressure and temperature, respectively, with ${\rm{\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 the temperature and pressure is found after the passing of the wavefront, i.e., some of the energy of the wave has been transformed into the internal energy of the plasma. The increase of the pressure seems to be uniform along the plasma. In contrast, it can be checked that the increase of the temperature is inversely proportional to the variation of the density.

The results shown in Figure 7 and in 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 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 density. In contrast, in partially ionized plasmas, the propagating waves leave a pressure wake due to the frictional dissipation of energy because of ion–neutral collisions, a phenomenon that is obviously absent from the fully ionized case of Verwichte et al. (1999).

Not all of the kinetic energy of the initial perturbation is used in heating the plasma, but a fraction of it is inverted in generating the second-order propagating waves. Hence, it is interesting to investigate how the energy deposition depends on the properties of the initial perturbation. A series of simulations has 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.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{e}_{k}(t=0)=\displaystyle \frac{1}{2l}{\int }_{-l}^{l}\displaystyle \sum _{s}{\rho }_{s}(x,t=0){[{V}_{y}^{(1)}(x,t=0)]}^{2}\,{dx},\end{eqnarray} \tag{ 5 }$$

and the variation of the internal energy of the medium is given by

$$\begin{eqnarray}&&{\rm{\Delta }}{e}_{P}(t)=\displaystyle \frac{1}{2l}{\int }_{-l}^{l}\displaystyle \sum _{s}\displaystyle \frac{{P}_{s}^{(2)}(x,t)-{P}_{s,0}(x)}{\gamma -1}\,{dx}.\end{eqnarray} \tag{ 6 }$$

The 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 ${\rm{\Delta }}{e}_{p}/{e}_{k}(t=0)\approx 6 \%$ is found at $\mathrm{FWHM}={10}^{5}\,{\rm{m}}$, which corresponds to a perturbation with an amplitude of ${V}_{y,0}=0.1{c}_{{\rm{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 with 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 the neutrals initially at rest. In another series of simulations, we perturb the y-component 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

$$\begin{eqnarray}&&{e}_{B}(t=0)=\displaystyle \frac{1}{2{l}_{1}}{\int }_{-{l}_{1}}^{{l}_{1}}\displaystyle \frac{{[{B}_{y}^{(1)}(x,t=0)]}^{2}}{2{\mu }_{0}}\,{dx}.\end{eqnarray} \tag{ 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 ∼2% when neutrals are initially at rest instead of ∼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: readers are reminded that under the chosen prominence conditions, neutrals account for two-thirds of the total mass of the plasma.

When the perturbation is applied to the y-component 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 FWHM ≈ 10⁵ m, and it has the same value as in the left panel of Figure 9, ∼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 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 the generation of waves and heating of the plasma. Hence, the internal energy has two components: one associated with the propagating wavefronts and another one related to the 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 energies is displayed for four simulations.

Zoom In Zoom Out Reset image size

In Figure 10, the total energy is not constant but decreases 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 ≈ 10 s, when most of the kinetic and magnetic energies goes to zero because Alfvén waves start crossing the boundaries. Later, the nonlinearly generated 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 ≈ 80 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 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 ≈ 10 s is larger in the bottom-right panel, which corresponds to an amplitude of ${V}_{y,0}=0.2{c}_{{\rm{A}}}$ and $\mathrm{FWHM}=1.25\times {10}^{4}\,{\rm{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 with 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.2{c}_{{\rm{A}}}$ than for ${V}_{y,0}=0.1{c}_{{\rm{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 nonlinear effects are more relevant. The reason is that more energy is required to generate second-order waves when the amplitude of the first-order perturbation increases, which leaves a smaller fraction of the initial energy that can be used to heat the plasma.

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 x-axis. 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 y-component of the velocity of all species for a time t_D, after which the driver is switched off. Hence, the driver can be cast as

$$\begin{eqnarray}{V}_{y,s}(x=0,t)=\left\{\begin{array}{lc}{V}_{y,0}\sin (\omega t) & \mathrm{if}\ t\leqslant {t}_{D},\\ 0 & \mathrm{if}\ t\gt {t}_{D,}\end{array}\right.\end{eqnarray} \tag{ 8 }$$

where the amplitude of the driver is ${V}_{y,0}=0.1{c}_{{\rm{A}}}$ and its frequency is $\omega ={10}^{-5}{{\rm{\Omega }}}_{p}\approx 1\,\mathrm{rad}\,{{\rm{s}}}^{-1}$, with ${{\rm{\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 over three 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 z-component, 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 y-components.

The figure displays the results after the driver has been switched off. The top-left panel shows the driven first-order wave, which propagates at the modified Alfvén speed, ${\widetilde{c}}_{{\rm{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 faster one travels with velocity ${\widetilde{c}}_{{\rm{A}}}$, and its longitudinal velocity is always positive. The slower 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\leqslant {\widetilde{c}}_{S}t$ results from the combination of the two waves. The frequency of the two waves is twice the frequency of the first-order 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\leqslant {\widetilde{c}}_{S}t$ while it increases in ${\widetilde{c}}_{S}\lt x\lt {\widetilde{c}}_{{\rm{A}}}t$. This means that the driver creates a mass flow that displaces the plasma away from the boundary. The aforementioned characteristics of the second-order perturbations are consistent with the analytical results shown in Appendix A.2, which have been derived for the simpler case of a two-fluid partially ionized plasma with strong ion–neutral 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 increase of the internal energy of the plasma due to the collisional friction. Although the driver used in this simulation fulfills $\omega /(2\pi )\lt {\nu }_{\mathrm{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\lt x\lt {\widetilde{c}}_{{\rm{A}}}t$ of the top-right panel, where ${V}_{x,s}(x,t)\gt 0$, indicating that there is mass flow toward 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 is 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.025{c}_{{\rm{A}}}$, ${V}_{y,0}=0.05{c}_{{\rm{A}}}$, and ${V}_{y,0}=0.1{c}_{{\rm{A}}}$, respectively. There is a slight decrease of the internal energy around t = 40 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.1{c}_{{\rm{A}}}$ is approximately four times larger than the relative variation for the case with ${V}_{y,0}=0.05{c}_{{\rm{A}}}$, and approximately 16 times larger than that for ${V}_{y,0}=0.025{c}_{{\rm{A}}}$, which corresponds to a quadratic dependence on the amplitude of the driver, as expected.

Zoom In Zoom Out Reset image size

The right panel shows the perturbations of the longitudinal component of the velocity of ions at a given time of the simulations. The plot focuses 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 that large if compared with ${\widetilde{c}}_{{\rm{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.1{c}_{{\rm{A}}}$ and ${V}_{y,0}=0.05{c}_{{\rm{A}}}$ are shown. As the leading perturbation travels to the right, the propagation speeds of the trailing waves change. One of the reasons for this change is the fluctuations in density. For instance, in the region ${\widetilde{c}}_{S}(t=0)t\lt x\lt {\widetilde{c}}_{{\rm{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\gt {\widetilde{c}}_{{\rm{A}}}(t=0)t$. On the other hand, ${\rm{\Delta }}\rho /{\rho }_{0}\lt 0$ for $x\lt {\widetilde{c}}_{S}(t=0)t$ and, consequently, ${\widetilde{c}}_{{\rm{A}}}(x,t)\gt {\widetilde{c}}_{{\rm{A}}}(t=0)$. The propagation speed of the sound waves is also modified by those density fluctuations, but it is also affected by the variation of the 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 first-order Alfvén wave due to ion–neutral collisions. Therefore, the second-order acoustic waves turn into shocks more easily in partially ionized plasmas than in fully ionized plasmas. Figure 13 shows that the fluctuations of ${\widetilde{c}}_{{\rm{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.

Zoom In Zoom Out Reset image size

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 ion–neutral 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 the entire running time of the simulation and not only for a brief time interval. Our goal here is to study the heating rate due to ion–neutral collisions. For this series of simulations, the driver is given by

$$\begin{eqnarray}&&{B}_{y,s}(x=0,t)={B}_{y,0}\sin (\omega t),\end{eqnarray} \tag{ 9 }$$

with ${B}_{y,0}={10}^{-2}{B}_{0}$. Thus, the energy input due to the effect of the driver is given by

$$\begin{eqnarray}&&{e}_{B}=\displaystyle \frac{{B}_{y}{(x=0,t)}^{2}}{2{\mu }_{0}}=\displaystyle \frac{{B}_{y,0}^{2}}{2{\mu }_{0}}{\sin }^{2}(\omega t),\end{eqnarray} \tag{ 10 }$$

and the mean energy input per period is

$$\begin{eqnarray}&&\langle {e}_{B}\rangle =\displaystyle \frac{1}{\tau }{\int }_{0}^{\tau }{e}_{B}\,{dt}=\displaystyle \frac{1}{4}\displaystyle \frac{{B}_{y,0}^{2}}{2{\mu }_{0}}.\end{eqnarray} \tag{ 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 through the perturbation of 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

$$\begin{eqnarray}&&{ \mathcal P }=\displaystyle \frac{\langle {e}_{B}\rangle }{\tau }=\displaystyle \frac{\omega \langle {e}_{B}\rangle }{2\pi },\end{eqnarray} \tag{ 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 larger power.

The heating rate can be computed as

$$\begin{eqnarray}&&Q=\displaystyle \sum _{s}\displaystyle \sum _{t\ne s}{Q}_{s}^{\mathrm{st}},\end{eqnarray} \tag{ 13 }$$

where Q_s^st is given by Equation (1) of Paper II. For the case of the three-fluid prominence plasma, the heating rate is then given by

$$\begin{eqnarray}Q & = & {\alpha }_{p{\rm{H}}}{({{\boldsymbol{V}}}_{{\boldsymbol{p}}}-{{\boldsymbol{V}}}_{{\rm{H}}})}^{2}+{\alpha }_{p\mathrm{He}}{({{\boldsymbol{V}}}_{{\boldsymbol{p}}}-{{\boldsymbol{V}}}_{\mathrm{He}})}^{2}\\ & & +{\alpha }_{\mathrm{HHe}}{({{\boldsymbol{V}}}_{{\rm{H}}}-{{\boldsymbol{V}}}_{\mathrm{He}})}^{2}+\displaystyle \frac{{\alpha }_{{pe}}}{{({{en}}_{e})}^{2}}{{\boldsymbol{j}}}^{2}\\ & & +{\alpha }_{e{\rm{H}}}{\left({{\boldsymbol{V}}}_{{\rm{H}}}-{{\boldsymbol{V}}}_{{\boldsymbol{p}}}+\displaystyle \frac{{\boldsymbol{j}}}{{{en}}_{e}}\right)}^{2}+{\alpha }_{e\mathrm{He}}{\left({{\boldsymbol{V}}}_{\mathrm{He}}-{{\boldsymbol{V}}}_{{\boldsymbol{p}}}+\displaystyle \frac{{\boldsymbol{j}}}{{{en}}_{e}}\right)}^{2},\end{eqnarray} \tag{ 14 }$$

where ${\boldsymbol{j}}={\sum }_{s}{Z}_{s}{{en}}_{s}{{\boldsymbol{V}}}_{{\boldsymbol{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}{{\rm{\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.

Zoom In Zoom Out Reset image size

The top panel of Figure 14 shows the normalized y-component 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 x-axis, 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 (-2{k}_{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 spatially averaged heating rate at every time step of the simulation in the following way:

$$\begin{eqnarray}&&{Q}_{\mathrm{av}}(t)=\displaystyle \frac{1}{{\lambda }_{I}}{\int }_{0}^{{\lambda }_{I}}Q(x,t)\,{dx}.\end{eqnarray} \tag{ 15 }$$

During the first time steps, Q_av grows as the waves propagate along the x-axis. 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.

Zoom In Zoom Out Reset image size

Song & Vasyliūnas (2011) found that the heating has a quadratic dependence on ω when $\omega \ll {\nu }_{\mathrm{st}}$ and that it is independent of ω 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 power-law function ${Q}_{\mathrm{fit}}\sim {\omega }^{a}$. In the range of low frequencies, the fitting exponent is a ≈ 1.93, quite close to the quadratic dependence found by Song & Vasyliūnas (2011). The exponent for intermediate frequencies switches to a ≈ 1.32. Finally, for higher oscillation frequencies, heating becomes independent of ω. 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, ω is much lower than all collision frequencies. Therefore, all of the species are strongly coupled and heating tends to a quadratic dependence on ω, as Song & Vasyliūnas (2011) found. However, as ω 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 two-fluid case of Song & Vasyliūnas (2011). Finally, at very high oscillation frequencies, heating does not depend on ω, in agreement with Song & Vasyliūnas (2011).

Figure 15 also allows the heating rate to be compared 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 ion–neutral collisions is more efficient when ω ≈ ν_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⁻⁵ to ${10}^{-4}\,{\rm{W}}\,{{\rm{m}}}^{-3}$. Comparing those values with the results represented in Figure 15, we find that L < Q for $\omega \lt {10}^{-4}\,{{\rm{\Omega }}}_{p}$. At low frequencies, heating by ion–neutral collisions can only compensate for a fraction of the radiated energy. In the range between $\omega \approx {10}^{-4}\,{{\rm{\Omega }}}_{p}$ and $\omega \approx {10}^{-3}\,{{\rm{\Omega }}}_{p}$, Q becomes of the order of L. At frequencies higher than $\omega \approx {10}^{-3}\,{{\rm{\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.

In this paper, we have continued the investigation started in Papers I and II about the importance of multi-fluid 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 this 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 multi-fluid model in which the effects of elastic collisions between all species of the plasma are taken into account. The general properties of nonlinear low-frequency Alfvén waves analyzed here are consistent with the results obtained by, e.g., Hollweg (1971), Rankin et al. (1994), Tikhonchuk et al. (1995), and Verwichte et al. (1999) for fully ionized plasmas, although differences appear due to the collisional interaction between ions and neutrals. For example, a second-order effect of nonlinear standing Alfvén waves is the appearance of a ponderomotive force that induces fluctuations in the density and pressure, and in 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 $2{k}_{z}{c}_{{\rm{A}}}$ and $2{k}_{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}}_{{\rm{A}}}$, and the weighted mean sound speed, ${\widetilde{c}}_{S}$, respectively, when the small-wavenumber range is considered. Since in the plasmas that have been examined here ${\widetilde{c}}_{{\rm{A}}}$ is much lower than ${\widetilde{c}}_{S}$, the second-order oscillations induced by the ponderomotive force are dominated by the mode associated with the sound speed. Due to this ponderomotive force, the plasma matter 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 multi-fluid 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 toward the anti-nodes 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 second-order 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) and Verwichte et al. (1999).

The propagation of nonlinear pulses through a plasma with conditions akin to those of a quiescent solar prominence has also been examined. The simulations have shown that after the initial perturbation has been applied to the plasma, the pulse splits into two smaller pulses that propagate in opposite directions at a speed given by ${\widetilde{c}}_{{\rm{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 the smaller width of the pulse. According to the results from Paper II, waves with larger wavenumbers have shorter damping times due to ion–neutral collisions, while perturbations with smaller wavenumbers are more long-lived. 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 second-order 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}}_{{\rm{A}}}$ and ${\widetilde{c}}_{S}$, respectively, as one of them is associated with 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 the form of heat. Consequently, the temperature of the plasma rises. Numerical simulations show that the heating generally increases when the width of the pulse is decreased, which, as already mentioned, is associated with efficient dissipation at small scales. However, at small enough widths, the calculated 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 second-order 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 first-order wave. For conditions in quiescent solar prominences, the investigation presented here has found that a maximum of 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 second-order perturbations in the density and pressure, and in the longitudinal component of the velocity, which propagate at the speeds ${\widetilde{c}}_{S}$ and ${\widetilde{c}}_{{\rm{A}}}$. As time advances, those perturbations cause the density to decrease in the region $x\leqslant {\widetilde{c}}_{S}t$ and to increase in the region ${\widetilde{c}}_{S}t\lt x\leqslant {\widetilde{c}}_{{\rm{A}}}t$. This means that the plasma is displaced toward 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 ion–neutral collisions enhance the formation of shocks from the second-order sound waves.

We have computed the heating rates due to elastic collisions and studied their spatial distribution and 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 (-2{k}_{I}x)$. On the other hand, when the oscillation frequency of the wave is much lower than the collision frequencies, the heating rates tend to a quadratic dependence on ω. 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 estimate of the radiative losses in a prominence ($L\sim {10}^{-5}\mbox{--}{10}^{-4}\,{\rm{W}}\,{{\rm{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 elaborate future studies that will exploit the full applicability of the nonlinear multi-fluid code developed here.

We acknowledge the support from grant AYA2014-54485-P (AEI/FEDER, UE). D.M. acknowledges support from MINECO through an FPI grant. R.S. acknowledges the Ministerio de Economía, Industria y Competitividad, and the Consellería d'Innovació, Recerca i Turisme del Govern Balear (Pla de ciència, tecnología, innovació i emprenedoria 2013-2017) for the "Ramón y Cajal" grant RYC-2014-14970. J.T. acknowledges support from MINECO and UIB through a "Ramón y Cajal" grant.

Here, a partially ionized two-fluid 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 ignored in 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

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{i}}{\partial t}+{\rm{\nabla }}\cdot ({\rho }_{i}{{\boldsymbol{V}}}_{{\boldsymbol{i}}})=0,\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{n}}{\partial t}+{\rm{\nabla }}\cdot ({\rho }_{n}{{\boldsymbol{V}}}_{{\boldsymbol{n}}})=0,\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial ({\rho }_{i}{{\boldsymbol{V}}}_{{\boldsymbol{i}}})}{\partial t} & + & {\rm{\nabla }}\cdot ({\rho }_{i}{{\boldsymbol{V}}}_{{\boldsymbol{i}}}{{\boldsymbol{V}}}_{{\boldsymbol{i}}})=-{\rm{\nabla }}{P}_{\mathrm{ie}}\\ & + & \displaystyle \frac{{\rm{\nabla }}\times {\boldsymbol{B}}}{{\mu }_{0}}\times {\boldsymbol{B}}+{\alpha }_{{in}}({{\boldsymbol{V}}}_{{\boldsymbol{n}}}-{{\boldsymbol{V}}}_{{\boldsymbol{i}}}),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial ({\rho }_{n}{{\boldsymbol{V}}}_{{\boldsymbol{n}}})}{\partial t}+{\rm{\nabla }}\cdot ({\rho }_{n}{{\boldsymbol{V}}}_{{\boldsymbol{n}}}{{\boldsymbol{V}}}_{{\boldsymbol{n}}})=-{\rm{\nabla }}{P}_{n}+{\alpha }_{{in}}({{\boldsymbol{V}}}_{{\boldsymbol{i}}}-{{\boldsymbol{V}}}_{{\boldsymbol{n}}}),\end{eqnarray} \tag{ 19 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\rm{\nabla }}\times ({{\boldsymbol{V}}}_{{\boldsymbol{i}}}\times {\boldsymbol{B}}),\end{eqnarray} \tag{ 20 }$$

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 ion–neutral friction coefficient. The expression for the friction coefficient can be found in Paper II. The rest of the symbols have already been defined.

To study the properties of nonlinear waves, a perturbative expansion is performed. Thus, each variable, ${\boldsymbol{f}}$, in the previous system of equations is rewritten as follows:

$$\begin{eqnarray}&&{\boldsymbol{f}}={{\boldsymbol{f}}}^{(0)}+\epsilon {{\boldsymbol{f}}}^{(1)}+{\epsilon }^{2}{{\boldsymbol{f}}}^{(2)}+...,\end{eqnarray} \tag{ 21 }$$

where 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 first-order and second-order perturbations, respectively. Since a static uniform background is considered, ${{\boldsymbol{V}}}_{{\boldsymbol{i}}}^{(0)}={{\boldsymbol{V}}}_{{\boldsymbol{n}}}^{(0)}=0$, and the remaining background values are constant.

Then, the terms in Equations (16)–(20) can be gathered according to their powers of , and separate 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 x-direction), and they are allowed to propagate along that same direction, the first-order (or linear) system leads to the equation for Alfvén waves,

$$\begin{eqnarray}&&\left[\displaystyle \frac{{\partial }^{3}}{\partial {t}^{3}}+(1+\chi ){\nu }_{{ni}}\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}-{c}_{{\rm{A}}}^{2}\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}-{c}_{{\rm{A}}}^{2}{\nu }_{{ni}}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}\right]{{\boldsymbol{V}}}_{{\boldsymbol{i}},\perp }^{(1)}=0,\end{eqnarray} \tag{ 22 }$$

where $\chi ={\rho }_{n}/{\rho }_{i}$ is the ionization ratio, ${\nu }_{{ni}}={\alpha }_{{in}}/{\rho }_{n}$ is the neutral–ion collision frequency, and ${{\boldsymbol{V}}}_{{\boldsymbol{i}},\perp }^{(1)}\equiv {V}_{i,y}^{(1)}\hat{{\boldsymbol{j}}}+{V}_{i,z}^{(1)}\hat{{\boldsymbol{k}}}$ is the perturbation of the velocity of ions in the direction perpendicular to the background magnetic field. The first-order perturbation of the magnetic field can be found through the equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{B}}}_{\perp }^{(1)}}{\partial t}={B}_{0}\displaystyle \frac{\partial {{\boldsymbol{V}}}_{{\boldsymbol{i}},\perp }^{(1)}}{\partial x},\end{eqnarray} \tag{ 23 }$$

where ${B}_{0}\equiv {B}_{x}^{(0)}$ is the background magnetic field and ${{\boldsymbol{B}}}_{\perp }^{(1)}\,\equiv {B}_{y}^{(1)}\hat{{\boldsymbol{j}}}+{B}_{z}^{(1)}\hat{{\boldsymbol{k}}}$.

The solutions of Equation (22) 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), and 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 the Alfvén and sound waves. In contrast, a coupling appears at the second order, as shown by the following equations, which are related to the velocities in the longitudinal direction:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{n}^{(2)}}{\partial t}+{\rho }_{n}^{(0)}\displaystyle \frac{\partial {V}_{n,x}^{(2)}}{\partial x}=0,\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{i}^{(2)}}{\partial t}+{\rho }_{i}^{(0)}\displaystyle \frac{\partial {V}_{i,x}^{(2)}}{\partial x}=0,\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{\rho }_{n}^{(0)}\displaystyle \frac{\partial {V}_{n,x}^{(2)}}{\partial t}=-\displaystyle \frac{\partial {P}_{n}^{(2)}}{\partial x}+{\alpha }_{{in}}({V}_{i,x}^{(2)}-{V}_{n,x}^{(2)}),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{\rho }_{i}^{(0)}\displaystyle \frac{\partial {V}_{i,x}^{(2)}}{\partial t}=-\displaystyle \frac{\partial {P}_{\mathrm{ie}}^{(2)}}{\partial x}-\displaystyle \frac{\partial }{\partial x}\left(\displaystyle \frac{{B}_{\perp }^{2}}{2{\mu }_{0}}\right)+{\alpha }_{{in}}({V}_{n,x}^{(2)}-{V}_{i,x}^{(2)}),\end{eqnarray} \tag{ 27 }$$

where ${B}_{\perp }^{2}\equiv {({B}_{y}^{(1)})}^{2}+{({B}_{z}^{(1)})}^{2}$. Thus, the second-order perturbation of the velocity of ions is related to the first-order 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 second-order equations corresponding to the perpendicular components have the same form as those of the first-order equations, and hence, they describe the same behavior as Equations (22) and (23).

The sound speeds of the ionized and neutral fluids are defined as ${c}_{\mathrm{ie}}=\sqrt{\gamma {P}_{\mathrm{ie}}^{(0)}/{\rho }_{i}^{(0)}}$ and ${c}_{S,n}=\sqrt{\gamma {P}_{n}^{(0)}/{\rho }_{n}^{(0)}}$, respectively. In the fully ionized single-fluid case, the second-order perturbations of pressure and density are related by the expression ${P}_{\mathrm{ie}}^{(2)}={c}_{\mathrm{ie}}^{2}{\rho }^{(2)}$ (see, e.g., Hollweg 1971; Rankin et al. 1994). When multi-fluid 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 (24)–(27), it is possible to obtain the following equation that describes the second-order perturbations of the density of ions (a similar equation can be cast for neutrals and for the x-component of the velocities of ions and neutrals):

$$\begin{eqnarray}&&\left[\displaystyle \frac{{\partial }^{4}}{\partial {t}^{4}}+({\nu }_{{in}}+{\nu }_{{ni}})\displaystyle \frac{{\partial }^{3}}{\partial {t}^{3}}-({c}_{S,n}^{2}+{c}_{\mathrm{ie}}^{2})\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}\right.\\ &&-\,\left.({\nu }_{{in}}{c}_{S,n}^{2}+{\nu }_{{ni}}{c}_{\mathrm{ie}}^{2})\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}+{c}_{\mathrm{ie}}^{2}{c}_{S,n}^{2}\displaystyle \frac{{\partial }^{4}}{\partial {x}^{4}}\right]{\rho }_{i}^{(2)}\\ &&=\,\left(\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}-{c}_{S,n}^{2}\displaystyle \frac{{\partial }^{4}}{\partial {x}^{4}}+{\nu }_{{ni}}\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}\right)\left(\displaystyle \frac{{B}_{\perp }^{2}}{2{\mu }_{0}}\right).\end{eqnarray} \tag{ 28 }$$

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,

$$\begin{eqnarray}&&\left[(1+\chi )\displaystyle \frac{{\partial }^{3}}{\partial {t}^{3}}-(\chi {c}_{S,n}^{2}+{c}_{\mathrm{ie}}^{2})\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}\right]{\rho }_{i}^{(2)}=\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}\left(\displaystyle \frac{{B}_{\perp }^{2}}{2{\mu }_{0}}\right),\end{eqnarray} \tag{ 29 }$$

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 on the right-hand side, namely

$$\begin{eqnarray}&&\left(\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}-{\widetilde{c}}_{S}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}\right){\rho }_{i}^{(2)}=\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}\left(\displaystyle \frac{{B}_{\perp }^{2}}{2{\mu }_{0}(1+\chi )}\right),\end{eqnarray} \tag{ 30 }$$

where an integration constant has been taken equal to zero, and ${\widetilde{c}}_{S}$ is an effective sound speed given by

$$\begin{eqnarray}&&{\widetilde{c}}_{S}={\left(\displaystyle \frac{{c}_{\mathrm{ie}}^{2}+\chi {c}_{S,n}^{2}}{1+\chi }\right)}^{1/2}.\end{eqnarray} \tag{ 31 }$$

Following a similar procedure, the differential equation for the second-order perturbation in the longitudinal component of the velocity is given by

$$\begin{eqnarray}&&\left(\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}-{\widetilde{c}}_{S}^{2}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}\right){V}_{i,x}^{(2)}=-\displaystyle \frac{1}{{\rho }_{i}}\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{\partial }{\partial x}\left(\displaystyle \frac{{B}_{\perp }^{2}}{2{\mu }_{0}(1+\chi )}\right).\end{eqnarray} \tag{ 32 }$$

From Equation (28), it is also possible to recover the differential equation that describes the second-order 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 (28) as

$$\begin{eqnarray}&&\left(\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}-{c}_{S,n}^{2}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}\right)\left(\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}-{c}_{\mathrm{ie}}^{2}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}\right){\rho }_{i}^{(2)}\\ &&\quad =\,\left(\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}-{c}_{S,n}^{2}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}\right)\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}\left(\displaystyle \frac{{B}_{\perp }^{2}}{2{\mu }_{0}}\right),\end{eqnarray} \tag{ 33 }$$

which leads to

$$\begin{eqnarray}&&\left(\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}-{c}_{\mathrm{ie}}^{2}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}\right){\rho }_{i}^{(2)}=\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}\left(\displaystyle \frac{{B}_{\perp }^{2}}{2{\mu }_{0}}\right),\end{eqnarray} \tag{ 34 }$$

an equation that has already been derived by Hollweg (1971), Tikhonchuk et al. (1995), and Terradas & Ofman (2004). It can be seen that Equations (30) and (34) represent the same type of behavior, with differences appearing in the velocity of the propagation of waves and the amplitude of the driving term. These are two effects caused by the ion–neutral interaction.

A.1. Standing Waves

In this section, the properties of nonlinear standing waves in a two-fluid partially ionized plasma are analyzed.

If the initial perturbation applied to the equilibrium state is given by

$$\begin{eqnarray}&&{V}_{y}^{(1)}(x,t=0)={V}_{y,0}\cos ({k}_{x}x),\end{eqnarray} \tag{ 35 }$$

and the strongly coupled limit is applied to Equations (22) and (23), the first-order perturbation of the magnetic field is

$$\begin{eqnarray}&&{B}_{\perp }(x,t)=\displaystyle \frac{-{B}_{0}}{{\widetilde{c}}_{{\rm{A}}}}{V}_{y,0}\sin ({\widetilde{c}}_{{\rm{A}}}{k}_{x}t)\sin ({k}_{x}x),\end{eqnarray} \tag{ 36 }$$

with ${\widetilde{c}}_{{\rm{A}}}$ the Alfvén speed modified by the inclusion of the inertia of neutrals, i.e., ${\widetilde{c}}_{{\rm{A}}}={B}_{0}/\sqrt{{\mu }_{0}{\rho }_{i,0}(1+\chi )}$. Then, using the initial conditions ${\rho }_{i}^{(2)}(x,t=0)=0$ and $\tfrac{\partial }{\partial t}{\rho }_{i}^{(2)}(x,t=0)=0$, respectively, the solution to Equation (30) can be computed as

$$\begin{eqnarray}&&{\rho }_{i}^{(2)}(x,t)=\displaystyle \frac{1}{2{\widetilde{c}}_{S}}{\int }_{0}^{t}{\int }_{x-{\widetilde{c}}_{S}(t-\tau )}^{x+{\widetilde{c}}_{S}(t-\tau )}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}\left[\displaystyle \frac{{B}_{\perp }^{2}(\xi ,\tau )}{2{\mu }_{0}(1+\chi )}\right]d\xi d\tau .\end{eqnarray} \tag{ 37 }$$

The only speed that explicitly appears in Equation (37) is the effective sound speed, ${\widetilde{c}}_{S}$. However, since the driving wave is assumed to be Alfvénic, ${B}_{\perp }^{2}$ is a function of the Alfvén speed. Hence, ${\rho }_{i}^{(2)}(x,t)$ also depends on that speed. Finally, the second-order perturbation of the density of ions is given by

$$\begin{eqnarray}&&{\rho }_{i}^{(2)}(x,t)\\ &&=\,\displaystyle \frac{{B}_{0}^{2}{V}_{y,0}^{2}[{\widetilde{c}}_{{\rm{A}}}^{2}-{\bar{c}}_{S}^{2}+{\bar{c}}_{S}^{2}\cos (2{\widetilde{c}}_{{\rm{A}}}{k}_{x}t)-{\widetilde{c}}_{{\rm{A}}}^{2}\cos (2{\bar{c}}_{S}{k}_{x}t)]\cos (2{k}_{x}x)}{8{\widetilde{c}}_{{\rm{A}}}^{2}{\bar{c}}_{S}^{2}({\widetilde{c}}_{{\rm{A}}}^{2}-{\bar{c}}_{S}^{2}){\mu }_{0}(1+\chi )}.\end{eqnarray} \tag{ 38 }$$

The resulting perturbation is the combination of two standing modes with frequencies $2{\widetilde{c}}_{{\rm{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, and ${\widetilde{c}}_{{\rm{A}}}$ with ${c}_{{\rm{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 (38) can be approximated as

$$\begin{eqnarray}&&{\rho }_{i}^{(2)}(x,t)\approx \displaystyle \frac{{B}_{0}^{2}{V}_{y,0}^{2}}{8{\mu }_{0}{\widetilde{c}}_{{\rm{A}}}^{2}(1+\chi )}\left[\displaystyle \frac{1-\cos (2{\widetilde{c}}_{S}{k}_{x}t)}{{\widetilde{c}}_{S}^{2}}\right]\cos (2{k}_{x}x),\end{eqnarray} \tag{ 39 }$$

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 1–3 is represented as ${\rm{\Delta }}\rho /{\rho }_{0}$, can be computed as the ratio between the second-order perturbation and the background density. Hence,

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}{\rho }_{i}}{{\rho }_{i,0}}\equiv \displaystyle \frac{{\rho }_{i}^{(2)}(x,t)}{{\rho }_{i,0}}\approx \displaystyle \frac{{V}_{y,0}^{2}}{8{\widetilde{c}}_{S}^{2}}[1-\cos (2{\widetilde{c}}_{S}{k}_{x}t)]\,\cos (2{k}_{x}x).\end{eqnarray} \tag{ 40 }$$

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}_{\mathrm{ie}}^{2}$ for fully ionized fluids and ${c}_{\mathrm{ie}}\gt {\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 that for smaller wavenumbers. Hence, it would be expected that Equation (40) 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 three-fluid simulations represented in Figures 2 and 3 show that, under the chosen physical parameters, the friction due to ion–neutral collisions is more efficient than the acoustic modes in attenuating the Alfvénic waves. For instance, it can be checked that in Figure 3 the first-order Alfvén wave has almost disappeared after t = 0.5 s, but the second-order perturbation in the x-component of the velocity lasts for a longer time. In a two-fluid plasma, the oscillation frequency and damping rate of the remaining second-order wave may be obtained from Equation (28) in the following way. Since the driving wave, i.e., the first-order Alfvén wave, vanishes due to collisions, after a given time, the term on the right-hand side of Equation (28) 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 second-order perturbation of the density of ions can be expressed as

$$\begin{eqnarray}&&{\rho }_{i}^{(2)}\sim \exp [i(-\omega t+\kappa x)],\end{eqnarray} \tag{ 41 }$$

where, in this case, the wavenumber is twice the wavenumber of the original driving wave, i.e., κ = 2k_x. This procedure leads to the following dispersion relation,

$$\begin{eqnarray}{\omega }^{4} & + & i({\nu }_{{ni}}+{\nu }_{{in}}){\omega }^{3}-{\kappa }^{2}({c}_{S,n}^{2}+{c}_{\mathrm{ie}}^{2}){\omega }^{2}\\ & & -i{\kappa }^{2}({\nu }_{{in}}{c}_{S,n}^{2}+{\nu }_{{ni}}{c}_{\mathrm{ie}}^{2})\omega +{c}_{\mathrm{ie}}^{2}{c}_{S,n}^{2}{\kappa }^{4}=0,\end{eqnarray} \tag{ 42 }$$

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 two-species fluid in which only the collisional interaction between ions and neutrals is taken into account and the influence of magnetic fields is ignored (see, e.g., Vranjes & Poedts 2010). It coincides with Equation (9) from Vranjes & Poedts (2010) if the factors proportional to the electron–neutral collision frequency of that formula are ignored, 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 speaking, the dispersion relation, Equation (42) should not be applicable because the driver is still working. This is a consequence of the damping due to ion–neutral 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}}_{{\rm{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 (38) and (39), if the Alfvén speed is much larger than the sound speed, the second-order oscillation is dominated by the acoustic mode. Hence, the results from Equation (42) 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 y-component of the velocity,

$$\begin{eqnarray}&&{V}_{y}^{(1)}(0,t)={V}_{y,0}\sin (\omega t),\end{eqnarray} \tag{ 43 }$$

where ${V}_{y,0}$ is proportional to the Alfvén speed, ${c}_{{\rm{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 x-axis, Equations (22) and (23) give the following first-order perturbations of the velocity and magnetic field:

$$\begin{eqnarray}{V}_{y}^{(1)}(x,t)=\left\{\begin{array}{lc}{V}_{y,0}\sin \left[\omega \left(t-\tfrac{x}{{\widetilde{c}}_{{\rm{A}}}}\right)\right], & \mathrm{if}\,x\leqslant {\widetilde{c}}_{{\rm{A}}}t,\\ 0, & \mathrm{if}\,x\gt {\widetilde{c}}_{{\rm{A}}}t,\end{array}\right.\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}{B}_{\perp }(x,t)=\left\{\begin{array}{lc}-\tfrac{{B}_{0}{V}_{y,0}}{{\widetilde{c}}_{{\rm{A}}}}\sin \left[\omega \left(t-\tfrac{x}{{\widetilde{c}}_{{\rm{A}}}}\right)\right], & \mathrm{if}\ x\leqslant {\widetilde{c}}_{{\rm{A}}}t,\\ 0, & \mathrm{if}\ x\gt {\widetilde{c}}_{{\rm{A}}}t.\end{array}\right.\end{eqnarray} \tag{ 45 }$$

The expression for the magnetic field perturbation can be inserted into Equation (32) to obtain the second-order 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 ({\delta }_{1}\omega t-{\delta }_{2}x)+{A}_{2}$ and inserting this expression into Equation (32). This step allows the values of A₁, δ₁, and δ₂ to be obtained. The value of A₂ is computed by imposing that ${V}_{\mathrm{NH}}^{(2)}(0,0)=0$. Then, the inhomogeneous solution is given by

$$\begin{eqnarray}&&{V}_{I}^{(2)}(x,t)=\displaystyle \frac{{V}_{y,0}^{2}{\widetilde{c}}_{{\rm{A}}}}{4({\widetilde{c}}_{{\rm{A}}}^{2}-{\widetilde{c}}_{S}^{2})}\left\{1-\cos \left[2\omega \left(t-\displaystyle \frac{x}{{\widetilde{c}}_{{\rm{A}}}}\right)\right]\right\}.\end{eqnarray} \tag{ 46 }$$

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 [{\delta }_{3}\omega (t-x/{\widetilde{c}}_{S})]+{A}_{4}$ and imposing that the complete solution fulfill the boundary condition ${V}_{i,x}^{(2)}(0,t)=0$, which means that there is no mass inflow. Hence,

$$\begin{eqnarray}&&{V}_{H}^{(2)}(x,t)=\displaystyle \frac{{V}_{y,0}^{2}{\widetilde{c}}_{{\rm{A}}}}{4({\widetilde{c}}_{{\rm{A}}}^{2}-{\widetilde{c}}_{S}^{2})}\left\{\cos \left[2\omega \left(t-\displaystyle \frac{x}{{\widetilde{c}}_{S}}\right)\right]-1\right\}.\end{eqnarray} \tag{ 47 }$$

Under the physical conditions analyzed in this work, Alfvén waves propagate faster than sound waves, as shown by Figure 11. Consequently, the second-order perturbation of the longitudinal component of the velocity is given by

$$\begin{eqnarray}{V}_{i,x}^{(2)}(x,t)=\left\{\begin{array}{lc}{V}_{H}^{(2)}(x,t)+{V}_{I}^{(2)}(x,t), & \mathrm{if}\ x\leqslant {\widetilde{c}}_{S}t,\\ {V}_{I}^{(2)}(x,t), & \mathrm{if}\ {\widetilde{c}}_{S}t\lt x\lt {\widetilde{c}}_{{\rm{A}}}t,\\ 0, & \mathrm{if}\ x\gt {\widetilde{c}}_{{\rm{A}}}t.\end{array}\right.\end{eqnarray} \tag{ 48 }$$

The second-order perturbation of density can be computed through the continuity equation, i.e., Equation (25), taking into account that ${\rho }_{i}^{(2)}(0,0)=0$. This leads to

$$\begin{eqnarray}{\rho }_{i}^{(2)}(x,t)=\left\{\begin{array}{lc}{\rho }_{H}^{(2)}(x,t)+{\rho }_{I}^{(2)}(x,t), & \mathrm{if}\ x\leqslant {\widetilde{c}}_{S}t,\\ {\rho }_{I}^{(2)}(x,t), & \mathrm{if}\ {\widetilde{c}}_{S}t\lt x\lt {\widetilde{c}}_{{\rm{A}}}t,\\ 0, & \mathrm{if}\ x\gt {\widetilde{c}}_{{\rm{A}}}t,\end{array}\right.\end{eqnarray} \tag{ 49 }$$

where

$$\begin{eqnarray}&&{\rho }_{I}^{(2)}(x,t)=\displaystyle \frac{{\rho }_{i,0}{V}_{y,0}^{2}}{4({\widetilde{c}}_{{\rm{A}}}^{2}-{\widetilde{c}}_{S}^{2})}\left\{1-\cos \left[2\omega \left(t-\displaystyle \frac{x}{{\widetilde{c}}_{{\rm{A}}}}\right)\right]\right\}\end{eqnarray} \tag{ 50 }$$

and

$$\begin{eqnarray}&&{\rho }_{H}^{(2)}(x,t)=\displaystyle \frac{{\rho }_{i,0}{V}_{y,0}^{2}{\widetilde{c}}_{{\rm{A}}}}{4{\widetilde{c}}_{S}({\widetilde{c}}_{{\rm{A}}}^{2}-{\widetilde{c}}_{S}^{2})}\left\{\cos \left[2\omega \left(t-\displaystyle \frac{x}{{\widetilde{c}}_{S}}\right)\right]-1\right\}.\end{eqnarray} \tag{ 51 }$$

The solutions for the fully ionized case are recovered by substituting ${\widetilde{c}}_{{\rm{A}}}$ and ${\widetilde{c}}_{S}$ with ${c}_{{\rm{A}}}$ and c_ie, respectively. In that case, Equations (46) and (47) 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 (46), (47), (50), and (51) show that the oscillation frequency of the second-order 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}\lt {\widetilde{c}}_{{\rm{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 than or equal to zero, while ${\rho }_{I}^{(2)}$ is larger than 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\lt {\widetilde{c}}_{S}t$ decreases.

The analytical expressions given by Equations (48) and (49) assume that the driver is always on. Nevertheless, they are good approximations for the perturbations of the longitudinal component of velocity and density during the initial seconds of the simulation represented in Figure 11.