© ESO 2022

1. Introduction

Alfvén waves are thought to play an important role in coronal heating, in solar wind acceleration, and in the development of turbulence in solar wind plasmas. We refer the reader to the works of Hasegawa & Uberoi (1982), Cross (1988), and Cramer (2001) for a broad introduction on the topic of Alfvén waves. Circularly polarised Alfvén waves are an exact solution to magnetohydrodynamic (MHD) equations even when the amplitude of the wave is large. The fluid approach has been used in the past in the description of Alfvén waves for reasons of mathematical tractability; however, fluid theory predicts a dramatic dependence of the occurrence of non-linear processes on the value of the plasma-β (e.g., Wong & Goldstein 1986). For this reason kinetic theory throughout the Vlasov equation provides the most general description of a plasma (see Spangler 1989, for applications to non-linear Alfvén waves).

The analysis of Alfvén waves over the last decades has been mostly focused on propagating waves because of their academic interest and also because of their practical applications to both laboratory and space plasmas. Hollweg (1971) showed that linearly polarised Alfvén waves propagating parallel to the magnetic field are able to drive density enhancements due to gradients in the wave magnetic-field pressure. In contrast, circularly polarised low-frequency Alfvén wave propagating along a constant magnetic field are known to be parametrically unstable (e.g., Derby 1978; Goldstein 1978), a process that has received significant attention in the scientific community. It can be understood as a wave-wave interaction mechanism in which a large-amplitude wave or pump couples non-linearly with smaller-amplitude fluctuations according to the laws of energy and momentum conservation. Three types of parametric instabilities are found: modulational, beat, and decay (e.g., Hollweg 1994). The decay instability is of particular interest in the solar wind since it might provide the backward-propagating Alfvén waves that are thought to be necessary for the development of MHD turbulence. The analysis of these instabilities has been mostly carried out in the single-fluid approach and dispersive effects arising from the cyclotron motion of single-ion species have been also included in the analysis. Moreover, kinetic effects have also been considered in the problem of parametric instability by Araneda (1998), Araneda et al. (2007, 2008). These authors, using linear Vlasov theory and hybrid simulations, found that kinetic effects such as Landau damping (LD) break the degeneracy of the mode-coupling solutions of fluid theory and are important even for very low proton plasma-β.

On the contrary, standing Alfvén waves have received much less attention in the literature. They can be formed at any place where trapped or reflected Alfvén waves propagate in opposite directions. They can be viewed as the superposition of two identical counter-propagating waves. The most clear examples of standing Alfvén waves are found in the magnetospheres of planets such as the Earth, Jupiter, and Mercury (see, e.g., Cummings et al. 1969; Manners et al. 2018; Khurana 1993). Other examples are related to standing Alfvénic transverse loop oscillations due to energy releases produced by flares, often reported in the solar atmosphere (e.g., Nakariakov et al. 1999; Schrijver et al. 2002; Aschwanden et al. 2002). A nearby perturbation is able to impact laterally coronal loops and excite Alfvén-like waves travelling along these closed magnetic flux tubes. These waves are almost fully reflected back at the solar photosphere where the density is several orders of magnitude higher than that in the solar corona (see the recent review of Nakariakov et al. 2021). Laboratory plasma is another example where the presence of an initial Alfvénic pump can generate standing Alfvén waves. For example, Danielson et al. (2002) studied the generation and damping of standing ion-acoustic waves in Penning-Malmberg electron traps.

Using a single-fluid approach Chian & Oliveira (1994, 1996), and Oliveira & Chian (1996) studied the effects of parametric instabilities on standing waves with applications to the magnetosphere. These authors developed a theory of MHD parametric instabilities driven by a left- or right-hand circularly polarised standing Alfvén wave in a low plasma-β under very specific initial conditions. Israelevich & Ofman (2004) and Israelevich & Ofman (2011) using both the MHD approach and hybrid simulations investigated numerically the generation of an electric field parallel to the magnetic field, due to non-linear induced ion-acoustic waves produced by the presence of a standing Alfvénic pump. In a different context, Mottez (2015) also explored the generation of a parallel electric field using PIC (particle in cell) simulations when two counter-propagating Alfvén waves are excited. However, the mechanism proposed by this author is not related to the ponderomotive force that it is usually operating under such conditions. A ponderomotive force (see the review of Lundin & Guglielmi 2006) generates density enhancements periodically in time and space which are normally called cavities and depend on gradients in density or electric fields. The periodicity depends on the value of gas pressure of the system and the amplitude of these enhancements is related quadratically with the amplitude of the initial transverse excitation (see also Rankin et al. 1994; Tikhonchuk et al. 1995). Interestingly, ponderomotive forces can lead to separation of species, although this topic is not addressed in the present work since the quasi-neutrality condition is imposed for the plasma composed of kinetic protons and the treatment of electrons is that of a fluid. Recently, Martínez-Gómez et al. (2018) used the multi-fluid approach to model high-frequency waves in partially ionised plasmas similar to those in quiescent solar prominences. These authors show that non-linear Alfvén waves generate density and pressure perturbations and that the friction due to collisions dissipates a fraction of the wave energy, which is transformed into heat that raises the temperature of the plasma.

The main goal of the present paper is to understand, using the multi-fluid approach and hybrid simulations, the generation of ion-acoustic modes driven non-linearly by standing finite amplitude transverse waves in 1D. We show both analytically and numerically that in our set-up parametric instabilities are secondary and that the ponderomotive force is mainly responsible for the dynamics of the system. Our analysis includes kinetic effects and addresses the relevant question of how LD of ion-acoustic waves due to wave-particle interactions is affected by the presence of the Alfvénic driver. We further propose some simple applications of the results presented in this paper.

2. Basic plasma model configuration

We assume a quasi-neutral electron-proton plasma. The background magnetic field, B₀, is constant and pointing in the x-direction. The equations are solved in this direction only, although the velocity and the electric and magnetic field components in the y- and z-directions are retained in the description (i.e., we concentrate on purely parallel wave propagation to the magnetic field). We adopt the multi-fluid approach and use the standard massless fluid-electron model. Nevertheless, protons and electrons are allowed to have different temperatures (T_p and T_e), and this leads to the presence of LD in the system, especially in the situation when T_p ≥ T_e.

Our interest is focused on standing waves, and therefore periodic boundary conditions are applied at the edge of the spatial domain, where nodes in the three components of the velocity are imposed at t = 0 and maintained during the time evolution. This mimics the effect of line-tying. The domain is defined in the range 0 ≤ x ≤ L, where L is the size of the system.

3. Preliminaries: Ion-acoustic standing waves

We start describing the results for standing ion-acoustic waves that are eigenmodes of the configuration and not forced oscillations due to a driver. Since the longitudinal magnetic field is assumed constant, with zero transverse magnetic fluctuations, it plays a passive role in this section. The results are well known in the case of propagating waves, but for standing waves there are some points that need to be explained in detail, the most relevant being that a standing ion-acoustic wave is attenuated by LD.

In the multi-fluid approach each species, denoted by the subindex s, satisfies the continuity equation and the momentum equation,

$$\begin{array}{cc}& \frac{\partial n_s}{\partial t}+\frac{\partial (n_s{u}_{\mathit{xs}})}{\partial x}=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\partial {\mathit{u}}_s}{\partial t}+u_{\mathit{xs}}\frac{\partial {\mathit{u}}_s}{\partial x}=\frac{q_s}{m_s}(\mathit{E}+\frac{{\mathit{u}}_s\times \mathit{B}}{c})-\frac{1}{m_s{n}_s}\frac{\partial p_s}{\partial x}{\widehat{\mathit{e}}}_x,\hfill \end{array}$$(2)

where collisions are neglected and the variables have their usual meaning. The mass density of each species is defined as ρ_s = m_sn_s. In the present case we only consider protons and electrons, but we assume that electrons are massless. This simplifies the problem and allows a direct comparison with the hybrid numerical simulations described later. The previous equations are completed with Maxwell equations, and the closure of the system of equation in fluids is often done with the adiabatic gas law. These equations lead to the MHD equations when the single-fluid approximation is used.

For a constant horizontal magnetic field the linearised longitudinal momentum equation reduces for the parallel velocity component to the following equation:

$$\begin{array}{c}\hfill m_s\frac{\partial \delta u_{\mathit{xs}}}{\partial t}=q_s\delta E_x-\frac{1}{n_{0s}}\frac{\partial \delta p_s}{\partial x}·\end{array}$$(3)

Here δu_xs, δE_x, and δp_s represent the parallel velocity, parallel electric field, and the pressure fluctuations, respectively, while n_0s is the equilibrium density. We eliminate the longitudinal component of the parallel electric field by assuming massless electrons; therefore, the previous equation when applied to electrons simplifies to (q_e = −e)

$$\begin{array}{c}\hfill \delta E_x=-\frac{1}{e{n}_0}\frac{\partial \delta p_{\mathrm{e}}}{\partial x},\end{array}$$(4)

where we neglect the term on the left-hand side of Eq. (3). Imposing quasi-neutrality means that n_0p = n_0e ≡ n₀, and for the same reason the proton and electron density fluctuations must be the same: δn_p = δn_e ≡ δn. Equation (4) is then introduced in Eq. (3) for protons. Using the standard adiabatic law for the perturbations we have

$$\begin{array}{c}\hfill \frac{\delta p_s}{p_{0s}}={\gamma }_s\frac{\delta n_s}{n_{0s}}·\end{array}$$(5)

These equations are complemented with the continuity equation which implies that δu_xp = δu_xe, and therefore hereafter we refer to this variable as δu_x.

The resulting system of equations is

$$\begin{array}{cc}& \frac{\partial \delta u_x}{\partial t}=-(c_{\mathrm{sp}}^2+c_{\mathrm{se}}^2)\frac{1}{n_0}\frac{\partial \delta n}{\partial x},\hfill \end{array}$$(6)

$$\begin{array}{cc}& \frac{\partial \delta n}{\partial t}=-n_0\frac{\partial \delta u_x}{\partial x}·\hfill \end{array}$$(7)

Here we use the definition of the sound speeds, $c_{\mathrm{se}}^2={\gamma }_{\mathrm{e}}\frac{k_{\mathrm{B}}}{m_{\mathrm{p}}}{T}_{\mathrm{e}}$ and $c_{\mathrm{sp}}^2={\gamma }_{\mathrm{p}}\frac{k_{\mathrm{B}}}{m_{\mathrm{p}}}{T}_{\mathrm{p}}$, being γ_e and γ_p the corresponding values of the polytropic index for electrons and protons, and T_e and T_p the equilibrium temperatures of the two species.

We rewrite Eq. (6) as a single equation for δu_x,

$$\begin{array}{c}\hfill \frac{{\partial }^2\delta u_x}{\partial t^2}-C_{\mathrm{s}}^2\frac{{\partial }^2\delta u_x}{\partial x^2}=0,\end{array}$$(8)

where we introduce $C_{\text{s}}^2=c_{\text{se}}^2+c_{\text{sp}}^2$. This equation is solved applying initial and boundary conditions. Since we are interested in standing waves in a finite domain these conditions are

$$\begin{array}{cc}& \delta u_x(x,t=0)=f(x),\hfill \end{array}$$(9)

$$\begin{array}{cc}& \frac{\partial \delta u_x}{\partial t}(x,t=0)=0,\hfill \end{array}$$(10)

$$\begin{array}{cc}& \delta u_x(0,t)=\delta u_x(L,t)=0,\hfill \end{array}$$(11)

where f(x) is the spatial profile of the initial perturbation. The solution to the wave equation together with the previous initial and boundary conditions is known to be

$$\begin{array}{c}\hfill \delta u_x(x,t)={\displaystyle \sum _{m=1}^{\infty }{c}_{m}cos\left(C_{\mathrm{s}}\frac{m\pi }{L}t\right)sin\left(\frac{m\pi }{L}x\right),}\end{array}$$(12)

where c_m are the Fourier coefficients given by

$$\begin{array}{c}\hfill c_m=\frac{2}{L}{\displaystyle {\int }_0^{L}f(x)sin\left(\frac{m\pi }{L}x\right)\mathrm{d}x,m=1,2,\dots }\end{array}$$(13)

Equation (12) satisfies the expected boundary conditions (nodes) at x = 0 and x = L in velocity. It represents a standing solution and can be viewed as the superposition of two counter-propagating waves. The standing solution has zero propagation speed since the two counter-propagating waves have the same phase velocity but of opposite sign.

The dispersion relation that we obtain is simply

$$\begin{array}{c}\hfill \omega =k{C}_{\mathrm{s}},\end{array}$$(14)

which corresponds to ion-acoustic waves, and

$$\begin{array}{c}\hfill k=\frac{m\pi }{L},m=1,2,\dots \end{array}$$(15)

is the discrete wave vector that matches the boundary conditions.

For the density perturbation we also obtain a homogeneous wave equation whose solution is calculated using Eqs. (7) and (12):

$$\begin{array}{c}\hfill \delta n(x,t)=-\frac{n_0}{C_{\mathrm{s}}}{\displaystyle \sum _{m=1}^{\infty }{c}_{m}sin\left(C_{\mathrm{s}}\frac{m\pi }{L}t\right)cos\left(\frac{m\pi }{L}x\right).}\end{array}$$(16)

Hence, density perturbations are 90° out of phase in time and space with respect to velocity. The density fluctuation, contrary to the parallel velocity, has antinodes at x = 0 and x = L. These results, although simple, are useful to understand the situation when the Alfvénic driver is present.

The expression for the electric field is obtained from Eq. (4) using Eq. (16) and reads

$$\begin{array}{c}\hfill \delta E_x(x,t)=-\frac{B_0}{c}\frac{1}{{\omega }_{\mathrm{cp}}}\frac{c_{\mathrm{se}}^2}{C_{\mathrm{s}}}{\displaystyle \sum _{m=1}^{\infty }\left(\frac{m\pi }{L}\right)c_{m}sin\left(C_{\mathrm{s}}\frac{m\pi }{L}t\right)sin\left(\frac{m\pi }{L}x\right),}\end{array}$$(17)

where the equilibrium magnetic field (B₀) and the cyclotron frequency for protons (ω_cp) appear for normalisation purposes only.

3.1. Kinetic effects on ion-acoustic standing waves

The aim of this section is to introduce kinetic effects by implementing a relatively simple method of Araneda (1998) and Araneda et al. (2007) to incorporate LD in the multi-fluid system of equations, previously described. The procedure considers the closure between pressure and density moments relationship by a fully kinetic treatment of the ion-acoustic wave dispersion analysis via the Vlasov equation, instead of using the traditional polytropic closure of the fluid equations. This avoids the inadequacies of fluid models in describing the resonant wave-particle interaction effects and the main result is that it leads to a polytropic index γ that is complex; see Appendix A for details about the procedure to include the kinetic effects in the standing problem. The most important conclusion of this analysis is that for the standing situation, the complex polytropic index γ is exactly the same as in the single-propagating case. This means that a standing wave will show exactly the same attenuation due to LD as the equivalent single-propagating wave, at least when there is no initial Alfvénic pump in the system. As we demonstrate below, the presence of the pump driver also does not change the attenuation by LD of the ion-acoustic wave.

Therefore, as shown by Araneda (1998) for the propagating case, to include the kinetic effects in the standing problem we just need to use the derived dispersion relation of the multi-fluid case and change γ of protons by Eq. (A.16) to obtain a new dispersion relation where ω is complex and takes into account the attenuation by LD. Hence,

$$\begin{array}{c}\hfill {\omega }^2=k^2(c_{\mathrm{se}}^2+c_{\mathrm{sp}}^2)=k^2[{\gamma }_{\mathrm{e}}\frac{k_{\mathrm{B}}}{m_{\mathrm{p}}}{T}_{\mathrm{e}}+2({\xi }_s^2-\frac{1}{Z^{\prime }(\xi )})\frac{k_{\mathrm{B}}}{m_{\mathrm{p}}}{T}_{\mathrm{p}}],\end{array}$$(18)

where $\xi =\frac{\omega }{k{\mathit{w}}_{0\mathrm{p}}}$, ${\mathit{w}}_{0\mathrm{p}}=\sqrt{2{\kappa }_{\mathrm{B}}{T}_{\mathrm{p}}/m_{\mathrm{p}}}$ is the thermal velocity of protons, and for simplicity we assume that the parallel drift is zero. It is convenient to define $c_{\mathrm{se}}={\mathit{v}}_{\mathrm{A}}\sqrt{{\gamma }_{\mathrm{e}}{\beta }_{\mathrm{e}}/2}$, and we have that T_p/T_e = β_p/β_e. From Eq. (18) is not difficult to write an explicit form for the dispersion relation in terms of the Z function:

$$\begin{array}{c}\hfill Z^{\prime }(\xi )=\frac{2{\beta }_{\mathrm{p}}}{{\gamma }_{\mathrm{e}}{\beta }_{\mathrm{e}}}·\end{array}$$(19)

This is exactly the same dispersion relation that is obtained by purely kinematic calculations for parallel-propagating ion-acoustic waves based on the Vlasov method (e.g. Swanson 2003). It is known that this dispersion relation has an approximate solution when the expansion for large arguments of the Z function is used. Under this approximation the real and imaginary parts of the frequency read

$$\begin{array}{cc}& {\omega }_{\mathrm{R}}=k{c}_{\mathrm{se}}\sqrt{1+3\frac{{\beta }_{\mathrm{p}}}{{\gamma }_{\mathrm{e}}{\beta }_{\mathrm{e}}}},\hfill \end{array}$$(20)

$$\begin{array}{cc}& {\omega }_{\mathrm{I}}=-k{c}_{\mathrm{se}}\sqrt{\frac{\pi {\gamma }_{\mathrm{e}}{\beta }_{\mathrm{e}}}{2{\beta }_{\mathrm{p}}}}{\xi }_{\mathrm{R}}^2{e}^{-{\xi }_{\mathrm{R}}^2},\hfill \end{array}$$(21)

where ${\xi }_{\mathrm{R}}=\frac{{\omega }_{\mathrm{R}}}{k{\mathit{w}}_{0\mathrm{p}}}$. This dimensionless parameter is found to be independent of the wavenumber and the sound speed using the approximation for ω_R given by Eq. (20). Interestingly, the damping per period defined as τ_D/P = ω_R/(2π|ω_I|) is also independent of k and c_se using the previous analytic approximations, and the dependence is only on the ratio β_p/β_e. It is worth mentioning that LD, apart from producing the attenuation of the oscillation, has a clear effect on the real part of the frequency of oscillation (see Eq. (20)), and this change in the period of oscillation might be easier to study in the time-dependent simulations rather than the damping time, which is more difficult to estimate based on the attenuation of the amplitudes with time. The implications of these results in the context of coronal loop oscillations are discussed in Sect. 5.

3.2. Hybrid simulations results

We performed hybrid simulations with the aim of understanding and corroborating the previous analytical results, especially regarding kinetic effects, and as a preparation for the driven case investigated below. The hybrid code we use is one-dimensional in space but retains all three components of the electromagnetic field, currents, and particle velocities (Winske & Leroy 1985; Terasawa et al. 1986; Horowitz et al. 1989; Winske & Omidi 1996). The electrons are considered as a massless and isothermal fluid, whereas ions are treated kinetically as discrete particles-in-cell (PICs) representing a collisionless plasma with finite generally anisotropic temperatures, and drifting parallel to the mean magnetic field. The advantage of using this hybrid method is that it allows us to resolve the ion dynamics and to integrate the equations over many ion-cyclotron periods, while neglecting the small temporal and spatial scales of the electron kinetic motions. The particle equations are advanced in time with a leapfrog method, and the fluid moments are computed with a second-order weighting scheme. The electric and magnetic fields are advanced in time solving the fields via a functional iterative explicit scheme, as described by Horowitz et al. (1989) and the derivatives in the gradient, curl, and divergence terms are estimated using a fourth-order finite difference scheme. A simple isothermal equation for electrons (γ_e = 1) is used. The code, modified from its original Fourier pseudo-spectral version by A. F. Viñas, has been tested and compared with similar numerical codes in several publications (e.g., Araneda et al. 2002, 2007; Xie et al. 2004; Moya et al. 2012; Maneva et al. 2013).

We use periodic boundary conditions in the system that mimic the effect of line-tying at the edges of the domain. We impose a sinusoidal initial perturbation in the parallel macroscopic velocity (i.e., a perturbation in δu_x) that excites a harmonic (k = 2 × 2π/L) of the ion-acoustic wave, which satisfies the periodic conditions. Hence, we are exciting a single mode in the system, which allows a comparison with the previous analytical results. Moreover, the chosen initial mode is precisely the mode that will be excited non-linearly in the driven case studied in Sect. 4. The macroscopic fluid velocity of the initial perturbation (u) is imposed self-consistently with the initial velocities of the individual particles (v) taking into account the shape of the macro-particles. Since we want to excite a standing mode, no density fluctuations need to be introduced in the system at t = 0 since density is phase-shifted with respect velocity in time by 90° in the standing wave, as we show above. This means that the initial position of the particles does not need to be imposed to obtain a certain initial density profile that excites an ion-acoustic wave. As usual, low-amplitude random velocities are introduced in the system. We typically use 256 grid cells with 1000 particles per cell. The time step is $\mathrm{\Delta }t=0.0025{\omega }_{\mathrm{cp}}^{-1}$. Lengths are normalised to the ion-inertial length, d_i = v_A/ω_cp and the velocity to the equilibrium Alfvén speed, v_A. We typically use a numerical box length of L = 20 d_i, although this parameter is not important regarding the damping per period of the ion-acoustic waves. We use that the ratio of the ion plasma frequency to the proton cyclotron frequency is ω_pi/ω_cp = 10⁴.

The time evolutions of different macroscopic variables are shown in Fig. 1 for a typical run. At the beginning of the simulation (top panel) the x-component of the current shows some coherent spatial structure, while the density fluctuation also displays a large scale associated with the mode number that is initially excited in the system (m = 4 in our case). We note that, as expected, J_x (essentially u_x) has nodes at the boundaries of the domain, while the density has antinodes. The magnetic field does not show a significant amplitude since the noise is low, and the same applies to the electric field components except for δE_x. As time evolves (bottom panel) the oscillatory character of the variables becomes visible and the standing character of the wave is clear. In these simulations we find a good agreement with the analytical expressions given previously. In particular, we checked that the overall structure of the parallel velocity and density is in agreement with Eqs. (12) and (16). Although the electric field seems quite noisy in Fig. 1, the longitudinal component behaves as predicted by Eq. (17). For a better comparison of the simulations and the multi-fluid solutions for the standing case we focused on a situation with very weak LD by choosing a small value of β_p in comparison with β_e.

Next we considered a situation with strong LD. The results of the simulations indicate that the amplitude of the excited ion-acoustic is attenuated with time. We inferred the corresponding period and damping time from the simulations. To do so we performed a Fourier analysis in space and focused on the dominant mode, which is precisely m = 4 because of the initial excitation. When Fourier analysis is performed the common definition of the wavenumbers is k = 2π/Lm′, and the mode m = 4 using Eq. (15) corresponds to the mode m′ = 2. An example of the temporal evolution associated with this dominant wavenumber is found in Fig. 2, top panel. The damping of the signal is quite clear in this plot, and the expected attenuation due to LD, by solving Eq. (19), is also plotted with a dashed envelope line for comparison purposes. The results from the dispersion relation are P = 44.6τ_A and τ_D = −127.2τ_A. From the simulations we performed a periodogram (Scargle 1982), which provides a better frequency accuracy than the common Fourier power spectrum, to calculate the real part of the frequency. An exponential fit to the envelope of the signal was used to estimate the imaginary part of the frequency. The obtained values from the simulations are P = 45.2τ_A and τ_D = −124.0τ_A. The agreement between the hybrid simulations and the theoretical result is rather good. Another example is shown in the bottom panel of Fig. 2. For the chosen parameters, with twice the plasma beta for protons, the attenuation is faster than in the previous case and the fit to the exponential decay is slightly worse. The results from the dispersion relation are P = 38.1τ_A and τ_D = −38.7τ_A, while the inferred values from the simulations are P = 36.9τ_A and τ_D = −42.9τ_A. For large time values the amplitude of the signal starts to be dominated by noise.

It is important to mention that the amplitude of the initial wave can have some influence on the attenuation of the signal. Since Landau theory is based on linear results, there might be some effects that are beyond the theoretical predictions. For this reason we decided to investigate the dependence of the attenuation with the amplitude of the initial wave using the hybrid simulations. The results are shown in Fig. 3 for two particular choices of values of (β_e, β_p). The variation in the real part of the frequency with the amplitude is quite small (top panel). While for the imaginary part the plot indicates that there is a sort of plateau where the inferred values tend to the theoretical calculations (plotted with horizontal dashed lines). Nevertheless, from Fig. 3 we realise that if the initial amplitude is too small or too large the deviations from the expected value can be significant. This is especially relevant for large values of the initial amplitude, and the reason behind the faster attenuation is that the excited ion-acoustic wave starts developing shocks. These shocks produce the excitation of higher spatial harmonics, and therefore some energy is deposited into these modes, meaning that the energy of the initial m = 4 mode decreases. Since the numerical code we are using does not contain proper shock-capturing techniques or explicit dissipation, this issue is not investigated further and we concentrate on the linear regime. On the other hand, for very small amplitudes, we have to remember that the system contains random noise and that in this case it is able to alter the efficiency of the LD process. We think that the results shown here are not related to the non-linear effects on the resonance itself. O’Neil (1965) showed that the damping can be significantly reduced because of the non-linear energy exchange between a wave and the resonant particles trapped in its potential wells. As the phase mixing of resonant particles becomes more complex, the energy interchange occurs less. In reality, it is known that the wave amplitude eventually saturates at a non-zero value which depends on the amplitude of the initial perturbation (e.g., Lancellotti & Dorning 1998, 2009).

Keeping in mind that the most convenient amplitude (typically around 0.01v_A) is in the proper attenuation regime dominated by linear LD, we constructed Fig. 4 by changing β_e and keeping the rest of the parameters constant. We represent again the real and imaginary part of omega together with the numerical solution of the dispersion relation and the analytical approximations. We find that the results inferred from the simulations agree quite well with the numerical solution for both the real and imaginary parts of ω. The analytical approximation is quite good for the real part of omega, while the imaginary part is less accurate (this is common in this kind of approximation; see e.g., Swanson 2003). The estimation of the damping time is more difficult from the simulations, especially when the attenuation is fast, but the fact that the real part agrees well with the theoretical prediction is already a clear indication that ion-acoustic waves lose their energy due to LD.

Further evidence of the presence of LD is found by calculating the velocity distribution function (VDF) from the hybrid simulations. In Fig. 5 we show the VDF as a function of the parallel velocity of the particles, v_x, and one of the transverse particle velocity components, v_y, integrated over the spatial simulation box along the x-direction. Similar results are obtained for the v_x − v_z velocity components. The initially circular profile at t = 0 (not shown in the plot) starts to develop two symmetrical beams along v_x since the standing wave can be viewed as the superposition of two identical counter-propagating waves. As time advances the beams are enhanced since more energy is deposited around the two symmetric resonance positions due to the wave-particle interaction. For a propagating wave we would only find a single beam instead of two, similar to the results obtained by Araneda et al. (2007).

The conclusion of this section is that standing ion-acoustic waves are eigenmodes of our finite length system and that these waves, similarly to the propagating case, attenuate by LD. In the linear regime the attenuation times and frequencies are exactly the same as those obtained for propagating waves and we can use the analytical approximations known for that situation. The fact that the standing wave is a combination of two counter-propagating waves does not change the efficiency of the wave-particle interaction process.

4. Ion-acoustic waves driven by Alfvén waves

The goal here is to investigate the generation of ion-acoustic waves driven non-linearly due to the presence of a standing Alfvénic pump. The previous section, based on the calculation of the non-forced ion-acoustic eigenmodes, turns out to be useful to understand the present driven case.

Since the Alfvénic standing pump can be viewed as the superposition of two circularly polarised counter-propagating waves, the first question that needs to be addressed is whether parametric instabilities are present in the system. For this reason we think that it is convenient to briefly re-examine the mathematical analysis that leads to the physical description of parametric instability for a circularly polarised propagating Alfvén wave. The idea is to use a perturbational approach and derive the equations at different orders in a dimensionless parameter. In this case (see Derby 1978; Goldstein 1978) the correct ordering in the perturbation scheme in ideal MHD equations is

$$\begin{array}{cc}& \mathit{u}={\mathit{u}}_{\perp }+ϵ{\mathit{u}}_{\perp }^{\prime }+ϵ\delta u_x{\widehat{\mathit{e}}}_x,\hfill \\ \hfill & \mathit{B}=B_0{\widehat{\mathit{e}}}_x+{\mathit{B}}_{\perp }+ϵ{\mathit{B}}_{\perp }^{\prime },\hfill \\ \hfill & \rho ={\rho }_0+ϵ\delta \rho ,\hfill \\ \hfill & p=p_0+ϵ\delta p\hfill \end{array}$$(22)

(the parameter is ϵ and is assumed to be much smaller than 1). We start with zero order in ϵ. We obtain from the perpendicular component of the momentum equation that

$$\begin{array}{c}\hfill {\rho }_0\frac{\partial {\mathit{u}}_{\perp }}{\partial t}=\frac{B_0}{{\mu }_0}\frac{\partial {\mathit{B}}_{\perp }}{\partial x},\end{array}$$(23)

while from the induction equation to zero order in ϵ we find

$$\begin{array}{c}\hfill \frac{\partial {\mathit{B}}_{\perp }}{\partial t}=B_0\frac{\partial {\mathit{u}}_{\perp }}{\partial x}·\end{array}$$(24)

Combining Eqs. (23) and (24) we obtain the standard wave equation that describes Alfvén waves. It is not difficult to see that if the propagating Alfvén wave is circularly polarised (the situation considered here), then B_⊥ = const; in other words, the amplitude of the wave is independent of time and position (see Appendix C) and this has important consequences, as we show in the following.

The parallel component of the momentum equation leads to zero order,

$$\begin{array}{c}\hfill 0=-\frac{\partial p_0}{\partial x}-\frac{1}{2{\mu }_0}\frac{\partial B_{\perp }^2}{\partial x},\end{array}$$(25)

which is a basic equation that is not clearly described in the literature of parametric instabilities. Since B_⊥ = const, this means that the second term in the previous equation (i.e. the ponderomotive term) is zero, and therefore p₀ must be constant as well and the equation is fulfilled. However, and this is the fundamental point of the present paper, for a standing circularly polarised wave the ponderomotive force of the pump is not zero and a different approach regarding the perturbation scheme needs to be used. For the moment we still provide the equations to first order in ϵ for the propagating circularly polarised wave. From the continuity equation we find

$$\begin{array}{c}\hfill \frac{\partial \delta \rho }{\partial t}+{\rho }_0\frac{\partial \delta u_x}{\partial x}=0.\end{array}$$(26)

The perpendicular component of the momentum equation gives

$$\begin{array}{c}\hfill {\rho }_0\frac{\partial {\mathit{u}}_{\perp }^{\prime }}{\partial t}+\delta \rho \frac{\partial {\mathit{u}}_{\perp }}{\partial t}+{\rho }_0\delta u_x\frac{\partial {\mathit{u}}_{\perp }}{\partial x}=\frac{B_0}{{\mu }_0}\frac{\partial {\mathit{B}}_{\perp }^{\prime }}{\partial x},\end{array}$$(27)

while for the parallel component

$$\begin{array}{c}\hfill {\rho }_0\frac{\partial \delta u_x}{\partial t}=-\frac{\partial \delta p}{\partial x}-\frac{1}{{\mu }_0}\frac{\partial ({\mathit{B}}_{\perp }·{\mathit{B}}_{\perp }^{\prime })}{\partial x}·\end{array}$$(28)

It is important to realise that on the right-hand side of this equation we find the coupling between the pump, B_⊥, and the induced magnetic perturbations, ${\mathit{B}}_{\perp }^{\prime }$. From the energy equation we obtain

$$\begin{array}{c}\hfill \frac{\partial \delta p}{\partial t}=-\gamma p_0\frac{\partial \delta u_x}{\partial x},\end{array}$$(29)

which is usually substituted in Eq. (28) to eliminate the variable δp. The magnetic field perturbation to first order in ϵ is

$$\begin{array}{c}\hfill \frac{\partial {\mathit{B}}_{\perp }^{\prime }}{\partial t}=B_0\frac{\partial {\mathit{u}}_{\perp }^{\prime }}{\partial x}-\frac{\partial (\delta u_x{\mathit{B}}_{\perp })}{\partial x},\end{array}$$(30)

where we see that there is a coupling with the longitudinal motions. These first-order equations in ϵ are exactly the same as those derived by Derby (1978) and Goldstein (1978) although this last author rewrites the equations as a second-order derivative in time and space. The same equations are found in Wong & Goldstein (1986) and Jayanti & Hollweg (1993) when dispersion effects, introduced by the finite ion cyclotron, are removed. These last authors, as in Derby (1978), prefer to write the equations in terms of B_± = B_y ± iB_z, which is useful when working with circularly polarised waves (the positive sign corresponds to left-handed waves and the negative sign to right-handed waves). It is well established that the previous equations lead to parametric instabilities (see the references in this paragraph).

If instead of a circularly polarised propagating wave we consider a standing circularly polarised wave, which is the main topic of the present work, we have that B_⊥ ≠ const (see Appendix C). This means that the perturbation expansion cannot be applied essentially because the ponderomotive force due to the pump is not zero now. The previous perturbation scheme fails precisely in Eq. (25) since the second term of the right-hand side is different from zero. Hence, we have to perform a different analysis to properly include the effect of the ponderomotive force in the equations, missing in the circularly polarised propagating case. It turns out that the appropriate expansion when B_⊥ ≠ const is the following (e.g., Rankin et al. 1994; Tikhonchuk et al. 1995; Ballester et al. 2020):

$$\begin{array}{cc}& \mathit{u}=ϵ{\mathit{u}}_{\perp }+{ϵ}^2\delta u_x{\widehat{\mathit{e}}}_x+{ϵ}^3{\mathit{u}}_{\perp }^{\prime },\hfill \\ \hfill & \mathit{B}=B_0{\widehat{\mathit{e}}}_x+ϵ{\mathit{B}}_{\perp }+{ϵ}^3{\mathit{B}}_{\perp }^{\prime },\hfill \\ \hfill & \rho ={\rho }_0+{ϵ}^2\delta \rho ,\hfill \\ \hfill & p=p_0+{ϵ}^2\delta p.\hfill \end{array}$$(31)

This expansion scheme is different from the situation for the circularly-propagating case given by Eq. (22). It is therefore logical to find a distinct physical behaviour in the two situations. In particular we realise that using the previous expansion the second-order density, pressure, and parallel velocity perturbations do not have the same order as the induced perpendicular velocity and magnetic field fluctuations which are third order.

We use the same procedure and collect terms of the same order in ϵ. From the perpendicular component of the momentum equation we obtain to first order

$$\begin{array}{c}\hfill {\rho }_0\frac{\partial {\mathit{u}}_{\perp }}{\partial t}=\frac{B_0}{{\mu }_0}\frac{\partial {\mathit{B}}_{\perp }}{\partial x},\end{array}$$(32)

while for the induction equation to first order in ϵ we find

$$\begin{array}{c}\hfill \frac{\partial {\mathit{B}}_{\perp }}{\partial t}=B_0\frac{\partial {\mathit{u}}_{\perp }}{\partial x}·\end{array}$$(33)

Combining Eqs. (32) and (33) we obtain again the standard wave equation that describes Alfvén waves. This is exactly the same result as in the previous propagating case, but now it is found to be to first order in ϵ. If we provide the equations to second order in ϵ, from the continuity equation we find

$$\begin{array}{c}\hfill \frac{\partial \delta \rho }{\partial t}+{\rho }_0\frac{\partial \delta u_x}{\partial x}=0.\end{array}$$(34)

From the parallel component of the momentum equation we obtain, to second order in ϵ that

$$\begin{array}{c}\hfill {\rho }_0\frac{\partial \delta u_x}{\partial t}=-\frac{\partial \delta p}{\partial x}-\frac{1}{2{\mu }_0}\frac{\partial B_{\perp }^2}{\partial x}·\end{array}$$(35)

This equation was first derived in Hollweg (1971; see his Eq. (11)) where a definite connection between linearly polarised Alfvén waves and density fluctuations was established. It is important to realise that on the right-hand side of this equation we encounter the effect of the ponderomotive force, and this equation is different from that obtained in the circularly polarised propagating case to first order (see second term on the right-hand side of Eq. (28)). This is what leads to different physical processes, namely a ponderomotive driven system in one case and to a parametrical unstable configuration in the other.

Interestingly, if we continue the analysis, to third order we obtain,

$$\begin{array}{c}\hfill {\rho }_0\frac{\partial {\mathit{u}}_{\perp }^{\prime }}{\partial t}+\delta \rho \frac{\partial {\mathit{u}}_{\perp }}{\partial t}+{\rho }_0\delta u_x\frac{\partial {\mathit{u}}_{\perp }}{\partial x}=\frac{B_0}{{\mu }_0}\frac{\partial {\mathit{B}}_{\perp }^{\prime }}{\partial x},\end{array}$$(36)

while for the magnetic field perturbation to third order in ϵ we have

$$\begin{array}{c}\hfill \frac{\partial {\mathit{B}}_{\perp }^{\prime }}{\partial t}=B_0\frac{\partial {\mathit{u}}_{\perp }^{\prime }}{\partial x}-\frac{\partial (\delta u_x{\mathit{B}}_{\perp })}{\partial x}·\end{array}$$(37)

These last two equations are exactly the same as in the circularly polarised propagating case (compare with Eqs. (27) and (30)). However, in the present scheme they are third order in ϵ, while the longitudinal velocity component is second order. This is different to the purely parametric unstable situation where density fluctuations and the induced velocity and magnetic field perturbations are of the same order in ϵ. Nevertheless, Eqs. (36) and (37) can still be viewed as the mathematical representation of a potential parametric process but to third order. Higher order terms can also play a role in regarding parametric instabilities. Thus, the main conclusion here is that for the standing wave problem parametric instabilities are not as primary as in the propagating case, and ponderomotive effects prevail. Therefore, we need to focus on the effect of the ponderomotive force, and for this reason we extend the previous results, using the perturbational approach given in Eq. (31) by using the multi-fluid equations that provide a more general description of a plasma (we assume again that electrons are massless). The electric field, present in the multi-fluid description, is assumed to have a similar expansion to the fluid velocity.

4.1. First-order linearisation

Collecting terms to first order in ϵ in the perturbational scheme we find the known equations that describe the circularly Alfvénic pump. Since we concentrate on standing waves instead of propagating waves we assume the following dependence (see Appendix C):

$$\begin{array}{cc}& u_{\perp }=i{u}_{\perp 0}{e}^{-i{\omega }_{0}t}sin{k}_{0}x,\hfill \end{array}$$(38)

$$\begin{array}{cc}& B_{\perp }=B_{\perp 0}{e}^{-i{\omega }_{0}t}cos{k}_{0}x.\hfill \end{array}$$(39)

We briefly describe here the features of this wave and how we obtain the dispersion relation and the polarisation relations in the multi-fluid case. The value of the wavenumber k₀ is discrete now.

Using Faraday’s law the electric field can be written in terms of the magnetic field as

$$\begin{array}{c}\hfill E_{\perp }=\frac{{\omega }_0}{c{k}_0}{B}_{\perp 0}{e}^{-i{\omega }_{0}t}sin{k}_{0}x.\end{array}$$(40)

The vectorial product between the velocity and the background is simply −iB₀u_⊥. To first order in ϵ we find

$$\begin{array}{c}\hfill -i{\omega }_0{u}_{\perp 0}=\frac{q_s}{m_s}(\frac{{\omega }_0}{c{k}_0}{B}_{\perp 0}+\frac{B_0}{c}\delta u_{\perp 0}).\end{array}$$(41)

If we introduce the cyclotron frequency ω_c = q_sB₀/m_sc (for protons q_s = e and m_s = m_p, while for electrons q_s = −e and m_s = m_e) the previous equation is written as

$$\begin{array}{c}\hfill ({\omega }_{\mathrm{c}}-{\omega }_0)u_{\perp 0}=-{\omega }_{\mathrm{c}}\frac{{\omega }_0}{k_0}\frac{B_{\perp 0}}{B_0}·\end{array}$$(42)

It is important to note that this equation is valid for both protons and electrons and that they have different velocity amplitudes in the perpendicular direction during the periodic oscillation. In the case of electrons (assuming that |ω_ce|≫ω₀) the relation reduces to

$$\begin{array}{c}\hfill u_{\perp 0\mathrm{e}}=-\frac{{\omega }_0}{k_0}\frac{B_{\perp 0}}{B_0}·\end{array}$$(43)

To derive a dispersion relation we have to use Ampere’s law which reads

$$\begin{array}{c}\hfill j_{\perp }=-i\frac{c}{4\pi }{k}_0{B}_{\perp 0}{e}^{-i{\omega }_{0}t}sin{k}_{0}x,\end{array}$$(44)

where we need the total perturbed current due to protons and electrons, which is

$$\begin{array}{c}\hfill j_{\perp }=e{n}_{0\mathrm{p}}{u}_{\perp \mathrm{p}}-e{n}_{0\mathrm{e}}{u}_{\perp \mathrm{e}}.\end{array}$$(45)

Due to charge neutrality the density of protons and electrons must be the same (i.e., n_0p = n_0e). Eliminating the velocities in favour of the magnetic field perturbation we have

$$\begin{array}{c}\hfill j_{\perp }=-e{n}_0\frac{B_{\perp }}{k_0{B}_0}\frac{{\omega }_0^2}{{\omega }_{\mathrm{cp}}-{\omega }_0}·\end{array}$$(46)

Substituting this expression in Eq. (44) we obtain the dispersion relation for standing waves which is exactly the same as that for propagating circularly polarised waves,

$$\begin{array}{c}\hfill 1-\frac{{\omega }_0}{{\omega }_{\mathrm{cp}}}=\frac{{\omega }_0^2}{{\omega }_{\mathrm{cp}}^2}\frac{{\omega }_{\mathrm{cp}}^2}{k_0^2{v}_{\mathrm{Ap}}^2},\end{array}$$(47)

where $v_{\text{Ap}}^2=B_0^2/(4\pi n_{\text{0p}}{m}_{\text{p}})$. Defining X₀ = ω₀/ω_cp and Y₀ = k₀v_Ap/ω_cp, the previous equations reduce to the standard form (e.g., Gomberoff & Elgueta 1991; Cramer 2001):

$$\begin{array}{c}\hfill Y_0^2=\frac{X_0^2}{1-X_0}·\end{array}$$(48)

For ω₀ > 0 this equation describes left-handed polarised Alfvén waves propagating in the positive direction (if k₀ > 0) or in the negative direction (if k₀ < 0) along the magnetic field, both with the same frequency.

For right-handed waves we need to change ω₀ by −ω₀ and the dispersion relation is

$$\begin{array}{c}\hfill Y_0^2=\frac{X_0^2}{1+X_0},\end{array}$$(49)

and again the same frequency is obtained for forward or backward waves since the dispersion relation is quadratic in k₀. Hence, we obtained a dispersion relation to first order for standing circularly polarised Alfvén waves which is exactly the same as for circularly polarised Alfvén propagating parallel to the magnetic field. The properties of the dispersion relation for L- and R-modes are shown in Fig. 6.

To finish this section it is interesting to realise that the polarisation relation given by Eq. (42) can be further simplified using the dispersion relation given by Eq. (47). We finally obtain that

$$\begin{array}{c}\hfill u_{\perp \mathrm{p}}=-\frac{v_{\mathrm{A}}^2}{{\omega }_0/k_0}\frac{B_{\perp }}{B_0}\end{array}$$(50)

for protons, while for electrons

$$\begin{array}{c}\hfill u_{\perp \mathrm{e}}=-\frac{{\omega }_0}{k_0}\frac{B_{\perp }}{B_0}·\end{array}$$(51)

Equations (50) and (51) are valid for both right- and left-handed circularly polarised modes, but ω₀ changes for R- and L-modes.

It should be noted that protons and electrons have the same velocity in MHD (i.e., when ω₀ = ±k₀v_A), and in this case

$$\begin{array}{c}\hfill u_{\perp }=-v_{\mathrm{A}}\frac{B_{\perp }}{B_0}\end{array}$$(52)

for k₀ > 0, and

$$\begin{array}{c}\hfill u_{\perp }=v_{\mathrm{A}}\frac{B_{\perp }}{B_0}\end{array}$$(53)

for k₀ < 0.

4.2. Second-order linearisation

To second order in ϵ we expect to obtain the dynamics of the wave due to ponderomotive effects. Using the standing Alfvénic pump to first order, described in the previous section, we have that for the parallel velocity component (same for protons and electrons) and to second order in ϵ

$$\begin{array}{c}\hfill m_s\frac{\partial \delta u_x}{\partial t}=q_s(\delta E_x+\frac{1}{c}{({\mathit{u}}_s\times \mathit{B})}_{2x})-\frac{1}{n_{0s}}\frac{\partial \delta p_s}{\partial x},\end{array}$$(54)

where to second order,

$$\begin{array}{c}\hfill {({\mathit{u}}_s\times \mathit{B})}_{2x}=\frac{1}{2i}(-u_{\perp s}{B}_{\perp }^{\ast }+u_{\perp s}^{\ast }{B}_{\perp }),\end{array}$$(55)

and protons and electrons have different perpendicular velocity amplitudes in the circularly polarised pump (see previous section for the relationship between u_⊥s and B_⊥).

As usual, we eliminate the longitudinal component of the parallel electric field by assuming massless electrons; therefore, from Eq. (54) applied to electrons and neglecting the term on the left-hand side, we obtain

$$\begin{array}{c}\hfill \delta E_x=-\frac{1}{c}{({\mathit{u}}_{\mathrm{e}}\times \mathit{B})}_{2x}-\frac{1}{e{n}_0}\frac{\partial \delta p_{\mathrm{e}}}{\partial x}·\end{array}$$(56)

This expression is then introduced in Eq. (54), but applied to protons.

Due to quasi-neutrality n_0p = n_0e = n₀, meaning also that δn_p = δn_e. Using the standard adiabatic law for the perturbations we have

$$\begin{array}{c}\hfill \frac{\delta p_s}{p_{0s}}=\gamma \frac{\delta n_s}{n_{0s}}·\end{array}$$(57)

We use this equation to eliminate the pressure perturbation and to work with the density perturbation.

The resulting system of equations using mass continuity reduces to

$$\begin{array}{cc}& \frac{\partial \delta u_x}{\partial t}=\frac{e}{m_{\mathrm{p}}}\frac{1}{c}(({\mathit{u}}_{\mathrm{p}}{\times \mathit{B})}_{2x}-{({\mathit{u}}_{\mathrm{e}}\times \mathit{B})}_{2x})-(c_{\mathrm{se}}^2+c_{\mathrm{sp}}^2)\frac{1}{n_0}\frac{\partial \delta n}{\partial x},\hfill \end{array}$$(58)

$$\begin{array}{cc}& \frac{\partial \delta n}{\partial t}=-n_0\frac{\partial \delta u_x}{\partial x}·\hfill \end{array}$$(59)

Here we use the definition of the sound speed, C_s.

We rewrite Eq. (58) as a single equation for δu_x,

$$\begin{array}{c}\hfill \frac{{\partial }^2\delta u_x}{\partial t^2}-(c_{\mathrm{se}}^2+c_{\mathrm{sp}}^2)\frac{{\partial }^2\delta u_x}{\partial x^2}=\frac{e}{m_{\mathrm{p}}}\frac{1}{c}\frac{\partial }{\partial t}({({\mathit{u}}_{\mathrm{p}}\times \mathit{B})}_{2x}-{({\mathit{u}}_{\mathrm{e}}\times \mathit{B})}_{2x}).\end{array}$$(60)

We see that we obtain a wave equation with a forcing term on the right-hand side due to the Alfvénic pump.

For the density we find the following equation:

$$\begin{array}{c}\hfill \frac{{\partial }^2\delta n}{\partial t^2}-(c_{\mathrm{se}}^2+c_{\mathrm{sp}}^2)\frac{{\partial }^2\delta n}{\partial x^2}=-\frac{n_0}{m_{\mathrm{p}}}\frac{e}{c}\frac{\partial }{\partial x}(({\mathit{u}}_{\mathrm{p}}{\times \mathit{B})}_{2x}-{({\mathit{u}}_{\mathrm{e}}\times \mathit{B})}_{2x}).\end{array}$$(61)

The term on the right-hand side now involves the spatial derivative of the forcing term instead of the time derivative that appears in the equation for the density fluctuation. Let us evaluate this forcing term. Using Eq. (38) for the standing pump it is not difficult to find that

$$\begin{array}{cc}& \frac{e}{m_{\mathrm{p}}}\frac{1}{c}(({\mathit{u}}_{\mathrm{p}}{\times \mathit{B})}_{2x}-{({\mathit{u}}_{\mathrm{e}}\times \mathit{B})}_{2x})\hfill \\ \hfill & =\frac{e}{m_{\mathrm{p}}}\frac{1}{c}(u_{\perp 0\mathrm{e}}-u_{\perp 0\mathrm{p}})B_{\perp 0}sin{k}_{0}xcos{k}_{0}x\hfill \\ \hfill & =\frac{{\omega }_{\mathrm{cp}}}{k_0}\frac{{\omega }_0^2}{{\omega }_{\mathrm{cp}}-{\omega }_0}\frac{B_{\perp 0}^2}{B_0^2}\frac{1}{2}sin2k_{0}x=Af(x),\hfill \end{array}$$(62)

where we write the velocity amplitudes for electrons and protons due to the pump wave in terms of the magnetic field perturbation. We rewrite the previous expression as a constant A times a function f(x). The forcing term is proportional to the square of the amplitude of the pump, as expected for a second-order expansion. However, we note that there is no temporal dependence. This means that the right-hand side of Eq. (60) is zero and the parallel motions seem to be uncoupled from the pump. Nevertheless, the forcing term is not zero for the density perturbation, and for this reason we concentrate on this magnitude. We solve Eq. (61) for the forcing term given by Eq. (62) together with the following initial and boundary conditions

$$\begin{array}{cc}& \delta n(x,t=0)=\frac{\partial \delta n}{\partial t}(x,t=0)=0,\hfill \end{array}$$(63)

$$\begin{array}{cc}& \frac{\partial \delta n}{\partial x}(0,t)=\frac{\partial \delta n}{\partial x}(L,t)=0.\hfill \end{array}$$(64)

It is known that the formal solution of this problem is the following,

$$\begin{array}{c}\hfill \delta n(x,t)={\displaystyle \sum _{m=1}^{\infty }[\frac{L}{C_{\mathrm{s}}m\pi }{\int }_0^t{f}_{m}sin\left(\frac{C_{\mathrm{s}}m\pi }{L}(t-\tau )\right)\mathrm{d}\tau ]cos\left(\frac{m\pi }{L}x\right),}\end{array}$$(65)

where

$$\begin{array}{c}\hfill f_m=\frac{2}{L}{\displaystyle {\int }_0^{L}g(x)cos\left(\frac{m\pi }{L}x\right)\mathrm{d}x,m=1,2,\dots }\end{array}$$(66)

We now have that

$$\begin{array}{c}\hfill g(x)=-n_{0}A\frac{\partial f(x)}{\partial x},\end{array}$$(67)

and we use that $C_{\mathrm{s}}=\sqrt{c_{\mathrm{se}}^2+c_{\mathrm{sp}}^2}$. Substituting f(x) into Eq. (67) and evaluating the integral in Eq. (66) we find, due to symmetry reasons, that all the values are zero except m = 4 (or m′ = 2) obtaining that

$$\begin{array}{c}\hfill f_4=-n_0{\omega }_{\mathrm{cp}}\frac{{\omega }_0^2}{{\omega }_{\mathrm{cp}}-{\omega }_0}\frac{B_{\perp 0}^2}{B_0^2}·\end{array}$$(68)

Using this expression and performing the corresponding integral with respect to τ in Eq. (65) we eventually obtain the following rather simple dependence for the perturbed density:

$$\begin{array}{c}\hfill \delta n(x,t)=n_0\frac{{\omega }_{\mathrm{cp}}}{{(C_{\mathrm{s}}2k_0)}^2}\frac{{\omega }_0^2}{{\omega }_{\mathrm{cp}}-{\omega }_0}\frac{B_{\perp 0}^2}{B_0^2}(cos\left(2k_0{C}_{\mathrm{s}}t\right)-1)cos(2k_{0}x).\end{array}$$(69)

This equation indicates that the characteristic wavenumber of the induced perturbations is 2k₀ and the corresponding frequency of the ion-acoustic driven wave is 2k₀C_s and it is independent of the frequency of the driver, ω₀. We can further simplify this expression using the dispersion relation given by Eq. (47) in the multiplicative factor, obtaining

$$\begin{array}{c}\hfill \delta n(x,t)=n_0\frac{v_{\mathrm{Ap}}^2}{4C_{\mathrm{s}}^2}\frac{B_{\perp 0}^2}{B_0^2}(cos\left(2k_0{C}_{\mathrm{s}}t\right)-1)cos(2k_{0}x).\end{array}$$(70)

Solving the wave equation for the velocity, using Eqs. (59) and (70) it is straightforward to find

$$\begin{array}{c}\hfill \delta u_x(x,t)=\frac{v_{\mathrm{Ap}}^2}{4C_{\mathrm{s}}}\frac{B_{\perp 0}^2}{B_0^2}sin\left(2k_0{C}_{\mathrm{s}}t\right)sin(2k_{0}x).\end{array}$$(71)

The previous expressions indicate a quadratic dependence of the density and velocity fluctuation with the perturbed perpendicular magnetic field. We recall that for a given perpendicular magnetic fluctuation associated with the Alfvénic wave the velocity amplitude for left- and right-handed modes is different according to, for example, Eq. (50).

It is worth mentioning that the previous result is for circularly polarised waves. For linearly polarised Alfvén waves an analytical solution is also available and it is different from the one obtained here, as expected. It was derived using the D’Alembert method (see Rankin et al. 1994; Tikhonchuk et al. 1995; Ballester et al. 2020). In this last case the parallel velocity is the sum of two sinusoidal terms with angular frequencies 2k₀c_s and 2 ω₀. Considering the low-β regime, c_s k₀ ≪ ω₀, then the dominant angular frequency of the solution is also 2k₀c_s.

It is useful to calculate an expression for the perturbed longitudinal electric field using Eq. (56) once we know the solution for the density perturbation, δn(x, t). We find that

$$\begin{array}{cc}\hfill \delta E_x(x,t)=& -\frac{1}{c}\frac{B_{\perp 0}^2}{B_0}\frac{1}{2}(\frac{{\omega }_0}{k_0}+\frac{c_{\mathrm{se}}^2}{C_{\mathrm{s}}^2}\frac{k_0}{{\omega }_{\mathrm{cp}}}{v}_{\mathrm{Ap}}^2)sin(2k_{0}x)\hfill \\ \hfill & +\frac{1}{c}\frac{B_{\perp 0}^2}{B_0}\frac{1}{2}\frac{c_{\mathrm{se}}^2}{C_{\mathrm{s}}^2}\frac{k_0}{{\omega }_{\mathrm{cp}}}{v}_{\mathrm{Ap}}^{2}cos\left(2k_0{C}_{\mathrm{s}}t\right)sin(2k_{0}x).\hfill \end{array}$$(72)

Hence, the parallel electric field has a time independent contribution (first term) plus another term that accounts for the ion-acoustic oscillation with time. At t = 0 the previous expression leads to

$$\begin{array}{c}\hfill \delta E_x(x,t=0)=-\frac{1}{c}\frac{B_{\perp 0}^2}{B_0}\frac{1}{2}\frac{{\omega }_0}{k_0}sin(2k_{0}x),\end{array}$$(73)

which is independent of the sound speed and only depends on the parameters of the Alfvén wave. At later stages of the evolution the generation of density fluctuations modifies the initial parallel electric field according to Eq. (72). However, we assume that there is a mechanism that attenuates the ion acoustic oscillation, for example, LD. If this is the case this type of oscillations will decay with time and its amplitude will tend to zero for t → ∞, but the Alfvénic wave will continue oscillating. We can use the previous expressions to find the ‘static state’ of the system. It is clear that if the ion-acoustic oscillation is attenuated with time then according to Eq. (71) it will eventually lead to δu_x(x, t) = 0. Although density is no longer uniform in the asymptotic state and using Eq. (70) it is found that

$$\begin{array}{c}\hfill \delta n(x,t\to \infty )=-n_0\frac{v_{\mathrm{Ap}}^2}{4C_{\mathrm{s}}^2}\frac{B_{\perp 0}^2}{B_0^2}cos(2k_{0}x).\end{array}$$(74)

The parallel electric field in this limit accounts only for the first term in Eq. (72):

$$\begin{array}{cc}\hfill \delta E_x(x,t\to \infty )& =-\frac{1}{c}\frac{B_{\perp 0}^2}{B_0}\frac{1}{2}(\frac{{\omega }_0}{k_0}+\frac{c_{\mathrm{se}}^2}{C_{\mathrm{s}}^2}\frac{k_0}{{\omega }_{\mathrm{cp}}}{v}_{\mathrm{Ap}}^2)sin\left(2k_{0}x\right)\hfill \\ \hfill & =-\frac{1}{c}\frac{B_{\perp 0}^2}{B_0}\frac{1}{2}\delta v_{\mathit{LR}}sin\left(2k_{0}x\right).\hfill \end{array}$$(75)

This means that there is a net parallel electric field in the system created by the Alfvénic pump that is able to sustain the created density cavities. These expressions provide information on how the background has changed once the induced ion-acoustic oscillations have disappeared due to some damping mechanism, for example LD. We define the velocity δv_LR as the term inside the brackets and this velocity is different for left- and right-handed modes. In Fig. 7 we see the dependence of δv_LR with the wavenumber for certain choice of values of plasma-β for electrons and protons. We conclude that the left-handed mode shows larger parallel electric field fluctuations (because v_LR is larger) than right-handed modes. The reason is that if ω₀ is positive for left-handed waves, but negative for right-handed waves then the addition of the two terms in the definition of v_LR (the second term is always positive) leading always to a larger value for the L-mode.

Hence, we solved analytically the problem of the non-linear excitation of the ion-acoustic due to standing circularly polarised Alfvénic pump in the absence of LD. Collisionless ion damping is addressed in a manner similar to that described in Sect. 3.1. The details of the formal derivation are found in Appendix D. The main result is that we can reproduce exactly the same procedure to introduce a kinetic effect in the multi-fluid systems via the polytropic relation between the pressure to second order and the density to second order. In other words, it is justified to substitute again the sound speed by its complex version that includes the effects of LD.

Now we apply the results of the second-order analysis to the standing case. Since we have seen that the pump simply excites the first harmonic with wavenumber 2k₀ of all the possible ion-acoustic eigenmodes, intuitively we expect that the frequency and damping times are precisely those found in Sect. 3.1 for 2k₀. However, we try to justify this result by assuming that ω is complex and substituting the polytropic γ_p for the corresponding complex expression involving the Z function. According to the analytical results in the absence of LD, we can assume that the velocity is of the form

$$\begin{array}{c}\hfill \delta u_x=Ai{e}^{-i\omega t}sin(2k_{0}x),\end{array}$$(76)

where ω can now be a complex number and A is a constant. From the continuity equation, given by Eq. (59) and using Eq. (76) we integrate with respect to time obtaining that

$$\begin{array}{c}\hfill \delta n=n_{0}Ai2k_{0}cos(2k_{0}x)\frac{1}{i\omega }(e^{-i\omega t}-1).\end{array}$$(77)

This expression is very similar to Eq. (70), but we have not determined the value of ω yet. To derive a dispersion relation we have to use the wave equation for the density given by Eq. (61) that includes the forcing term, but if we want to include the LD effects we have already shown that we have to use the complex version of C_s (complex γ_p), defined in Eq. (18). From the density wave equation we find, after cancelling factors that contain the spatial dependence, the following simple equation:

$$\begin{array}{c}\hfill \frac{1}{i\omega }2k_{0}Ai{(-i\omega )}^2{e}^{-i\omega t}+C_{\mathrm{s}}^2\frac{1}{i\omega }{(2k_0)}^{3}Ai(e^{-i\omega t}-1)=-k_0^2{v}_{\mathrm{Ap}}^2\frac{B_{\perp 0}^2}{B_0^2}·\end{array}$$(78)

This equation has two oscillatory terms proportional to e^−iωt and two non-oscillatory terms. To make the sum of the two oscillatory terms zero, because the equation must be valid for any time, we find the following condition must be satisfied,

$$\begin{array}{c}\hfill \omega =2k_0{C}_{\mathrm{s}},\end{array}$$(79)

where C_s depends on ω now and contains the Z function. This equation is precisely the dispersion relation for ion-acoustic waves that undergo LD. The corresponding wavenumber is twice that of the Alfvénic driver.

In addition, the non-oscillatory term with time on the left-hand side of Eq. (78) must be equal to the term on the right-hand side, meaning that

$$\begin{array}{c}\hfill A=\frac{v_{\mathrm{Ap}}^2}{4C_{\mathrm{s}}}\frac{B_{\perp 0}^2}{B_0^2},\end{array}$$(80)

where we use Eq. (79) to eliminate ω from the equation. Since we are working with complex numbers we have to take the real part of Eqs. (76) and (77). If we assume that C_s is real (and therefore independent of ω) we recover the results of Sect. 4.2, to be compared with the amplitude in Eq. (71).

The main conclusion is that we obtain exactly the same dispersion relation as in Sect. 3.1 for non-driven ion-acoustic waves. In the driven case we simply have to change k by 2k₀ and solve the dispersion relation or use the analytical expressions for the real and imaginary parts of the frequency already derived. The Alfvénic pump, at second order in ϵ, does not modify the attenuation by LD of the generated ion-acoustic waves. The reason is that the forcing term on the right-hand side of Eq. (78) does not contain any dependence with the sound speed, C_s, that accounts for LD effects.

4.3. Hybrid simulations results for the Alfvénic pump driver

As we described in Sect. 3.2, we performed hybrid simulations to have further confirmation of the previous analytical results based on the previous linear analysis. An example of a typical run is shown in Fig. 8. The choice of the ratio of the electron to proton temperatures makes the LD mechanism inefficient in this case and therefore suitable to be compared with the obtained analytical results. We chose to excite the system with a right-handed circularly polarised Alfvén wave with a wavelength that matches the size of our spatial domain; therefore, periodic boundary conditions are applicable (i.e. k₀ = 2π/L). In order to excite this type of wave we need to impose a certain profile for the velocity and magnetic field components at t = 0 (see details in Appendix C; in particular, we have to use the excitation given by Eqs. (C.13)–(C.16)). In the hybrid simulation the particle distribution function is initiated self-consistently with the bulk velocity of the pump driver. The velocity and the magnetic field perturbations are related (see Eq. (42)) through the fluid dispersion relation for the particular Alfvén driver pump wave (given by Eq. (47)). There is no initial fluid speed in the parallel direction and no initial density perturbation either (since the Alfvénic wave is incompressible). In our case the initial distribution contains a delta function with drifts representing the transverse part and multiplying a Gaussian function with no longitudinal drift so that

$$\begin{array}{c}\hfill f=\delta (v_{\perp }-u_{\perp 0})e^{-{(v_x-u_{x0})}^2/{w_0}^2}\end{array}$$(81)

with u_x0 = 0 (see also Eqs. (A.3) and (A.4)). This initial excitation is compatible, for example, with that of Sonnerup & Su (1967).

The simulation shows that we have excited the right-handed wave since its frequency agrees with the expected value, given by Eq. (49). We clearly find the driven ion-acoustic wave that is non-linearly excited in the system, and we are also able to reproduce its behaviour using the previous analytical results. For example, the density shows periodic spatial changes (cavities) with a wavenumber that is twice that of the Alfvénic driver (2k₀). We overplotted the analytical profiles for the parallel velocity and density fluctuation (Eqs. (70) and (71)) in Fig. 8 (see thick orange dashed lines) and we find that the agreement with the simulation results is quite remarkable. In addition, the parallel electric field closely matches the formula given by Eq. (72). Similar results are obtained if instead of exciting the right-handed wave we perturb the system with a left-handed polarised wave, see Eqs. (C.29)–(C.32).

We now analyse how the driver affects the LD of the excited wave. For comparison purposes we consider the same parameters in Sect. 3.2 where we studied the behaviour of an initially excited ion-acoustic wave. In the present simulations we find that the induced ion-acoustic shows an exponential damping, similar to that shown in Fig. 2. Nevertheless, we find some dependence of the attenuation on the amplitude of the Alfvénic pump and therefore on the amplitude of the ion-acoustic wave generated non-linearly. We already found this effect in the undriven case (see Sect. 3.2). We performed several simulations by changing the initial amplitude of the Alfvénic pump and constructed Fig. 9. We find that the smaller the amplitude of the pump, the weaker the attenuation by LD. Only when the initial velocity amplitude is around 0.1v_A does the damping per period approach the theoretical result of the non-driven ion-acoustic attenuated wave. However, we have to be cautious about this value of the pump since as we keep increasing its value more energy goes into the ion-acoustic wave (with an amplitude that is approximately given by Eq. (80)) which can eventually start developing shocks, as we found in Sect. 3.2. In any case, the agreement of the real part of omega with the expected value from the theoretical prediction is a clear indication that the LD is taking place and that the process of attenuation is not significantly altered in the driven case. This agrees with the analytical justification that we give in Sect. 4.2.

Finally, we investigate the generation of the parallel electric field in the simulations. We calculate the Fourier power associated to the m′ = 2 mode (k = 2k₀), which is the dominant wavenumber in the parallel electric field. The results from the hybrid calculations are shown in Fig. 10. We find that the electric field starts at a certain value, but then decreases and oscillates around a new value with the frequency of the driven ion-acoustic wave. This oscillation shows attenuation due to the LD process. The initial value of the parallel electric field agrees with the calculation given by Eq. (73) valid at t = 0 (see dotted line). The asymptotic value also matches the theoretical prediction for the stationary state given by Eq. (75) (see dashed line). This corroborates the validity of our theoretical analysis.

Finally, from the hybrid simulations, we investigated the possible presence of parametric instabilities. For most of the cases investigated in the present paper we do not find clear evidence of the presence parametric unstable modes. Only in some specific situations with very small LD damping and very long simulation times can the effects of parametric instabilities be identified. This confirms that the role of parametric instabilities is secondary in the present standing problem, as is demonstrated in the analysis given at the beginning of Sect. 4.

5. Application of the results

We briefly describe a possible application of the results discussed so far. Slow modes, closely related to ion-acoustic modes, are often reported in the solar corona. In particular, observations indicate that standing wave slow modes are commonly excited in coronal loops (see the review of Wang et al. 2021). These modes are identified by Doppler and intensity fluctuations in hot coronal lines. A possible mechanism that can explain the presence of these modes, at least the first harmonic, is precisely a transverse Alfvénic oscillation that eventually drives non-linearly ion-acoustic waves. One of the most interesting features of the reported oscillations is that they are heavily damped with time; the damping per period, τ_D/P, is of the order of one, meaning a severe attenuation. Several mechanisms, including non-adiabatic effects such as thermal conduction, compressive viscosity, and optically thin radiation have been proposed. Non-linearity, the cooling loop background, the wave-caused heating–cooling imbalance, plasma non-uniformities, and other effects such as loop geometry, wave leakage through footpoints, and corona and magnetic effects (mode coupling and obliqueness) have also been used to try to explain the fast attenuation. Voitenko et al. (2005) concluded that the observed dissipation distances for propagating slow modes can alternatively be explained by phase mixing in its ideal regime, where the apparent damping is due to the spatial integration of the phase mixed amplitudes by the observation.

The measurements of the difference in electron and ion temperature in the solar corona indicate in some cases that T_e ≈ T_p/3.5 (e.g., Boldyrev et al. 2020). This value suggests that ion-acoustic waves suffer strong LD under coronal conditions and therefore this wave-particle interaction mechanism may explain, for example, the reported attenuation of standing slow modes in coronal loops. The damping per period due to LD is essentially independent of the wavenumber (see Eqs. (20) and (21)), meaning that although we focus in the present work on very small spatial scales, in reality the mechanism operates on much larger spatial scales as well. Similar results concerning kinetic damping of small-scale slow-wave modes and higher frequency waves were obtained by Viñas et al. (2000) with regard to coronal heating and the formation of suprathermal tails in the coronal plasma. Therefore, we propose here that the attenuation of standing slow modes reported in the solar corona might be linked to the process of LD. The expression for the damping per period using Eqs. (20) and (21) is

$$\begin{array}{c}\hfill {\tau }_{\mathrm{D}}/P=\frac{\sqrt{2}}{{\pi }^{3/2}}{e}^{\frac{1}{2}(\psi +3)}\frac{1}{\psi }\frac{1}{\sqrt{\psi +3}},\end{array}$$(82)

where

$$\begin{array}{c}\hfill \psi ={\gamma }_{\mathrm{e}}\frac{{\beta }_{\mathrm{e}}}{{\beta }_{\mathrm{p}}}·\end{array}$$(83)

This expression is only applicable in the situation β_p ≤ β_e (i.e., for weak attenuation). Equation (82) can be used in the future for seismological purposes. In particular, it might be possible to infer the ratio β_e/β_p based on the reported characteristics of the damped ion-acoustic modes.

The effect of LD has been investigated, in the coronal context, for propagating waves (see Ruan et al. 2016a,b). These authors performed simulations of driven slow modes along the magnetic field by using a kinetic model that includes the effects of LD and Coulomb collisions. They found that the obtained spatial scale of damping is similar to that resulting from thermal conduction in MHD. Future studies should address the case of the standing wave problem studied here, but including the effect of collisions that are not negligible in the lower solar corona and that might have a significant effect on the attenuation by LD.

6. Conclusions and discussion

Two cases of the evolution of ion-acoustic waves have been investigated. We first studied the excitation of purely ion-acoustic standing modes using a multi-fluid approach. We solved analytically the initial and boundary value problem incorporating kinetic effects by using a complex polytropic index. We showed that ion-acoustic standing waves attenuate by LD, a result that is completely equivalent to the Landau attenuation for propagating ion-acoustic waves extensively explored in the literature. The performed numerical hybrid simulations have shown however that the attenuation follows the predictions of the LD up to a certain degree. Too large or too small initial excitation amplitudes do not lead to the expected damping times due to the effects of non-linearity and/or noise. In general, we find that for the standing case two symmetric beams develop in the distribution function since the standing wave can be viewed as the superposition of the two counter-propagating waves.

The presence of two symmetric counter-streaming beams in the distribution function can have important consequences regarding the stability of the system. For example, it is well known in plasma physics that counter-streaming beams can be the source of the oscillating two-stream instability that strongly affects the stability of such distributions and becomes the source for additional wave mode excitation. In the future this aspect should be investigated in detail since two stream instabilities may develop and affect the dynamics of the configuration. In essence, when two streams move through each other one wavelength in one cycle of the plasma frequency a density perturbation on one stream is reinforced by the forces due to bunching of particles in the other stream and vice versa (e.g., Dawson 1962; Birdsall & Langdon 1991).

The second aspect investigated in this paper considers the case when the Alfvénic pump is present, driving the ion-acoustic waves. For this case we showed that standing circularly polarised waves have the same dispersion relation as propagating waves, and we derived the appropriate polarisation conditions for their excitation. We further demonstrated that parametric instabilities are secondary in the standing problem, contrary to the situation studied by Chian & Oliveira (1994, 1996), Oliveira & Chian (1996), and to the purely propagating case. The Alfvénic pump, via ponderomotive forces, excites non-linearly the second harmonic of the standing ion-acoustic wave (i.e. with k = 2k₀, where k₀ is the wavenumber of the standing pump). The frequency of the ion-acoustic wave is ω = 2k₀C_s, in the absence of Landau attenuation. When LD is relevant it also affects the real part of the frequency. These ion-acoustic waves are damped with time depending on the ratio T_e/T_p. The hybrid simulations have shown that the attenuation follows the standard linear predictions in the non-driven case up to a certain degree. The attenuation is weaker when the amplitude of the driver is low.

We derived analytical expressions for the density and electric fluctuations induced quadratically by the pump. Left- and right-handed modes produce exactly the same density fluctuations, but we have demonstrated that the L-mode has the largest parallel electric field fluctuation. These results provide a possible interpretation of measurements of electric fields parallel to the magnetic field observed, for example, in the auroral ionosphere. Although several mechanisms can explain this reported electric field, Israelevich & Ofman (2004, 2011) have proposed that precisely standing Alfvén waves can naturally generate the parallel electric field. The analytic expressions for the generated electric fields could be used in the future to compare with the observations, although the first impression is that this electric field is rather small in general. On a different context, Mottez (2015) has investigated the generation of electric fields due to the interaction of counter-propagating Alfvén waves. Nevertheless, the generated electric field in our case is purely due to the ponderomotive forces, while the electric field described by the previous author has a completely different origin.

We propose that standing slow modes reported in coronal loops can be interpreted in terms of ion-acoustic waves that attenuate by LD. The ratio T_e/T_p plays the relevant role in this problem, and perhaps can be indirectly inferred from the damping per period of slow modes reported from the observations. However, the low corona is known to be weakly collisional, which can be a competing effect that can have a significant influence on the process of LD. The efficiency of the attenuation is altered by the thermalisation of the non-thermal tail in the velocity distribution by Coulomb collisions. In the propagating case beam flows are produced by Landau resonance, but they are destroyed by Coulomb friction. From the theoretical point of view the inclusion of collisional effects is through the Vlasov equation by including a Fokker-Planck collisional term, although there are other alternatives. This problem should be investigated in the near future also with potential applications to partially ionised plasmas.

The reasons to neglect electron LD are that the focus of this paper is to study low-frequency waves (i.e., waves below the ion-cyclotron and ion-plasma frequencies) for which electron LD has little or no contribution since this effect is more important for high-frequency waves (e.g., Whistler, electron-acoustic, and Langmuir waves). Furthermore, because the ion and electron temperatures are similar, the electron thermal velocity (i.e., the width) of the electron velocity distribution function is about 43 times larger than the ion thermal velocity. This indicates that the electron LD is much smaller than that of ion LD for resonances near the ion-thermal speed since LD is proportional to the slope of the velocity distribution function. Only ion LD is included in this work. In addition, the use of the hybrid simulation neglects by construction electron LD since electrons are treated as a fluid. Additionally, there are various reasons to consider only zero ion perpendicular temperatures in our study. This is due to the fact that a finite ion-perpendicular temperature has no effect on LD and only contributes to ion-cyclotron damping when finite ion-gyroradius effects are included. In principle, finite perpendicular temperature effects can be included; however, the mathematical calculations become more tediously complex since in this case one needs to introduce not only a single cyclotron frequency, but also an infinite series of cyclotron-harmonics that will require a truncation. More importantly, the effect of finite perpendicular temperature will not contribute to LD at all, only to cyclotron damping. Since in our study we are only considering parallel propagation, this further eliminates the effect of finite ion-perpendicular temperatures. The hybrid simulations include finite-temperature effects such as cyclotron-damping, but they do not change the LD, which is a longitudinal effect.

The results presented here will be used to consider more realistic situations. In particular, it is interesting to introduce different ion populations to study their effects. For example, two counter-propagating light beams together with a heavy stationary population. This could represent a loop with weak flows penetrating thought the footpoints. The stability of this system can be then analysed using kinetic theory with both analytical methods and hybrid simulations. Interestingly, different types of instabilities can appear in the system. Along the same line, we focused on initial excitations that correspond to pure standing waves, meaning that the two counter-propagating waves have exactly the same amplitude and wavelength but propagate in opposite directions. It would be interesting to analyse the effect of two slightly unbalanced counter-propagating waves (i.e., when the symmetry is broken). This can significantly change the dynamics of the system, and can lead to the appearance of clearer beams and possibly to the development of stronger parametric instabilities, which are weak in the perfectly symmetric case.

Movies

Movie 1 associated with Fig. 1 (fig1) Access here

Movie 2 associated with Fig. 8 (fig8) Access here

Acknowledgments

J.T. acknowledges the support from grant AYA2017-85465-P (MINECO/AEI/FEDER, UE), to the Conselleria d’Innovació, Recerca i Turisme del Govern Balear, and also to IAC³. A.F. Viñas would like to thank the Universitat de les Ilhes Balears for their assistance as a visiting professor and to the Catholic University/IACS and NASA-GSFC for their support during the development of this work. J.A. acknowledges the support from FONDECYT grant 1161700. The authors thank Prof. Marcel Goossens for useful suggestions that helped to improve the paper. We thank the anonymous referee for their useful comments that have helped to improve this manuscript.

References

Appendix A: Incorporating kinetic effects in the fluid description

The method was presented Araneda (1998) and Araneda et al. (2007) for propagating waves, but here we extend it to the case of standing waves. Kinetic effects in the fluid equations emerge through the Vlasov equation,

$$\begin{array}{cc}\hfill \frac{\partial F_s(x,\mathbf{v},t)}{\partial t}& +\mathbf{v}·\mathrm{\nabla }{F}_s(x,\mathbf{v},t)\hfill \\ \hfill & +\frac{q_s}{m_s}(\mathbf{E}(x,t)+\frac{\mathbf{v}\times \mathbf{B}(x,t)}{c})\frac{\partial F_s(x,\mathbf{v},t)}{\partial \mathbf{v}}=0,\hfill \\ \hfill \end{array}$$(A.1)

where F_s(x, v, t) represents the distribution function of species s and v is the velocity vector that corresponds to the particle distribution velocities. It is important not to confuse this velocity with the fluid velocity, u, introduced in the main text.

We linearise Eq. (A.1) with

$$\begin{array}{cc}\hfill F_s(x,\mathbf{v},t)& =F_{0s}(x,\mathbf{v},t)+\delta f_s(x,\mathbf{v},t),\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill F_{0s}(x,\mathbf{v},t)& =\delta ({\mathbf{v}}_{\perp }-{{\mathbf{u}}_{\perp }}_{0s}(x,t))f_{0s}(v_x),\hfill \end{array}$$(A.3)

where we assume that the fluctuation of the velocity distribution function is δf_s(x, v_x, t). In Eq. (A.3) we assume that the perpendicular temperature in the zeroth order is zero as represented by the Dirac delta-function in the transverse velocity direction, which includes only the effect of the transverse wave velocity drift. The zeroth-order parallel velocity distribution f_0s(v_x) can be taken, with all generality, as a finite temperature drifting Maxwellian distribution given by

$$\begin{array}{c}\hfill f_{0s}(v_x)=\frac{n_{0s}}{\sqrt{\pi }{w}_{0s}}{e}^{-\frac{{(v_x-u_{x0s})}^2}{w_{0s}^2}},\end{array}$$(A.4)

where n_0s is the equilibrium ion density, ${\mathit{w}}_{0s}=\sqrt{2{\kappa }_B{T}_{0s}/m_s}$ is the thermal velocity, and u_x0s is the velocity drift. The first-order linearised expression for δf_s(x, v_x, t) can be immediately determined as

$$\begin{array}{cc}\hfill \frac{\partial \delta f_s(x,v_x,t)}{\partial t}& +v_x\frac{\partial \delta f_s(x,v_x,t)}{\partial x}+\frac{q_s}{m_s}(\delta E_x(x,t)+\hfill \\ \hfill & \frac{\delta u_{\mathit{ys}}{B}_z-\delta u_{\mathit{zs}}{B}_y+u_{\mathit{ys}}\delta B_z-u_{\mathit{zs}}\delta B_y}{c})\frac{\partial f_{0s}(v_x)}{\partial v_x}=0.\hfill \\ \hfill \end{array}$$(A.5)

We assume that the fluctuations vary as e^{i(kx − ωt)} with a complex frequency ω and a real wave vector k. We therefore consider propagating waves for the moment, but as mentioned earlier the goal is to find a solution for standing waves. Following this linearising procedure in Eq. (A.5) yields

$$\begin{array}{c}\hfill \delta f_s(x,v_x,t)=-i\frac{q_s}{m_s}\delta G_s(x,t)\frac{1}{\omega -k{v}_x}\frac{\partial f_{0s}(v_x)}{\partial v_x},\end{array}$$(A.6)

where we define

$$\begin{array}{c}\hfill \delta G_s(x,t)=\delta E_x+\frac{\delta u_{\mathit{ys}}{B}_z-\delta u_{\mathit{zs}}{B}_y+u_{\mathit{ys}}\delta B_z-u_{\mathit{zs}}\delta B_y}{c},\end{array}$$(A.7)

which represents a longitudinal force; it is dependent on (x, t), but independent of the particle velocities. At this point is unnecessary to know the exact form of this force, thus we focus on the resonant wave-particle interaction term in Eq. (A.6). We note that for purely ion-acoustic waves propagating with parallel propagation, the magnetic terms are zero in Eq. (A.7) and only the x-component of the electric field survives.

We now calculate the moment density fluctuation δn_s of the plasma from the linearised velocity distribution fluctuation δf_s(x, v_x, t) as

$$\begin{array}{c}\hfill \delta n_s(x,t)={\displaystyle {\int }_{-\infty }^{+\infty }\delta f_s(x,v_x,t)d{v}_x.}\end{array}$$(A.8)

By substitution of Eq. (A.4) and Eq. (A.6) into Eq. (A.8) and performing the integrals (see the details in Appendix B) we obtain

$$\begin{array}{c}\hfill \delta n_s(x,t)=i{n}_{0s}\frac{q_s}{m_s}\frac{\delta G_s(x,t)}{k{w}_{0s}^2}{Z}^{\prime }({\xi }_s),\end{array}$$(A.9)

which is written in terms of the derivative of the Zeta function, Z, and ${\xi }_s=\frac{\omega -k{u}_{x0s}}{k{\mathit{w}}_{0s}}$. We now associate the density fluctuation derived in Eq. (A.9) with the pressure fluctuation via the usual polytrope equation

$$\begin{array}{c}\hfill \frac{\delta p_s}{p_{0s}}=\gamma \frac{\delta n_s}{n_{0s}},\end{array}$$(A.10)

where γ is the polytrope coefficient that is not constant anymore, but depends on the complex frequency, the wave vector, and the background equilibrium parameters of the plasma. To obtain the proper representation of the polytrope γ we need to use the linearised continuity equation and the linearised longitudinal fluid momentum equation for species s given by

$$\begin{array}{cc}\hfill \frac{\partial \delta n_s}{\partial t}& +\frac{\partial (\delta n_s{u}_{0xs}+n_{0s}\delta u_{\mathit{xs}})}{\partial x}=0,\hfill \end{array}$$(A.11)

$$\begin{array}{cc}\hfill \frac{\partial \delta u_{\mathit{xs}}}{\partial t}+& u_{0xs}\frac{\partial \delta u_{\mathit{xs}}}{\partial x}=\hfill \\ \hfill & \frac{q_s}{m_s}(\delta E_x+\frac{\delta u_{\mathit{ys}}{B}_z-\delta u_{\mathit{zs}}{B}_y+u_{\mathit{ys}}\delta B_z-u_{\mathit{zs}}\delta B_y}{c})\hfill \\ \hfill & -\frac{1}{m_s{n}_{0s}}\frac{\partial \delta p_s}{\partial x}.\hfill \end{array}$$(A.12)

Upon linearisation of Eq. (A.11) and Eq. (A.12) we obtain

$$\begin{array}{cc}\hfill (\omega -k{u}_{0xs})\delta n_s& =k{n}_{0s}\delta u_{\mathit{xs}},\hfill \end{array}$$(A.13)

$$\begin{array}{cc}\hfill -i(\omega -k{u}_{0xs})\delta u_{\mathit{xs}}& =\frac{q_s}{m_s}\delta G_s(x,t)-i\frac{k\delta p_s}{m_s{n}_{0s}}.\hfill \end{array}$$(A.14)

From Eq. (A.13) we get an expression for δu_xs, which is substituted in Eq. (A.14). Combined with Eq. (A.9) it leads to the complex kinetic polytropic relationship

$$\begin{array}{c}\hfill \frac{\delta p_s}{p_{0s}}=2({\xi }_s^2-\frac{1}{Z^{\prime }({\xi }_s)})\frac{\delta n_s}{n_{0s}},\end{array}$$(A.15)

where Z′ is the derivative of the plasma dispersion function. We also assumed that p_0s = n_0sk_BT_0s.

Therefore, comparing Eq. (A.10) with Eq. (A.15) we recognise the complex kinetic polytrope coefficient as

$$\begin{array}{c}\hfill \gamma =2({\xi }_s^2-\frac{1}{Z^{\prime }({\xi }_s)}).\end{array}$$(A.16)

Equations (A.15) and (A.16) clearly show that density and pressure fluctuations are not necessarily in phase, in contrast with the phase-locking that results from the constant value obtained in the fluid treatment. This result indicates that the kinetic polytropic coefficient depends on the dispersive properties of the plasma. Since this coefficient is in the Fourier-domain, its inverse transform (e.g. via a superposition of the wave modes that corresponds to multiple spatial scales of the system) of Eq. (A.16) implies an intricate non-local real space relationship between pressure and density of the system. On the contrary, in the fluid treatment the polytrope is assumed to be constant and the pressure and density relationship to be local. Therefore, the dispersive property of the polytrope coefficient introduces non-local effects which ultimately control the transfer of energy between finite amplitude transverse and longitudinal fluctuations (Mjolhus & Wyller 1988; Spangler 1989; Hammett & Perkins 1990). As a result this can explain, to a certain extent, why it is difficult to obtain reasonable values of the polytropic coefficient when comparing linear fluid theory with kinetic simulations (Vasquez 1995). We should note, however, that we included kinetics effects only in the parallel direction, but we kept the assumption for drifting Maxwellian distributions. This means that we keep only the collisionless damping effects of LD on the growth rates of ion-acoustic waves, but all other kinetic cyclotron resonance effects (due to the simplified zeroth transverse temperature effects and finite gyro-radius effects) are completely neglected.

We repeat the same analysis, but for a propagating wave in the opposite direction; the wave now propagates to the left, meaning that we are assuming that the wavenumber is −k. Following the same procedure we now have that the fluctuation in the perturbed distribution function is simply

$$\begin{array}{c}\hfill \delta f_s(x,v_x,t)=-i\frac{q_s}{m_s}\delta G_s(x,t)\frac{1}{\omega +k{v}_x}\frac{\partial f_{0s}(v_x)}{\partial v_x},\end{array}$$(A.17)

but we note the change in the denominator with respect to the equivalent equation for the right-propagating wave in Eq. (A.6).

Calculating the density fluctuation δn_s of the plasma from the linearised velocity distribution fluctuation δf_s(x, v_x, t) as before, we obtain now

$$\begin{array}{c}\hfill \delta n_s(x,t)=-i{n}_{0s}\frac{q_s}{m_s}\frac{\delta G_s(x,t)}{k{w}_{0s}^2}{Z}^{\prime }({\xi }_s),\end{array}$$(A.18)

where we use the definition for ξ_s as in the previous case (${\xi }_s=\frac{\omega }{k{\mathit{w}}_{0s}}$, where for simplicity we eliminated the drift). Details about this derivation are found in Appendix B. We note the change in the global sign in comparison with the result for the wave propagating to the right (see Eq. (A.9)). The calculation of the complex kinetic polytropic relationship now leads to

$$\begin{array}{c}\hfill \frac{\delta p_s}{p_{0s}}=2({\xi }_s^2-\frac{1}{Z^{\prime }({\xi }_s)})\frac{\delta n_s}{n_{0s}},\end{array}$$(A.19)

and therefore the complex kinetic polytrope coefficient is

$$\begin{array}{c}\hfill \gamma =2({\xi }_s^2-\frac{1}{Z^{\prime }({\xi }_s)}).\end{array}$$(A.20)

It is clear that from the physical point of view it should lead to exactly the same damping times as for the right-propagating wave when there is no Alfvénic pump. This agrees with the value of γ that we obtain; it is exactly the same as in Eq. (A.16) for the wave propagating to the right.

These results are for propagating waves, but we are interested in standing waves. The previous derivation leads us to think that in fact the full standing solution for the distribution function should be of the form

$$\begin{array}{c}\hfill \delta f_s(x,v_x,t)=f_{+}(v_x)e^{i(kx-\omega t)}+f_{-}(v_x)e^{i(-kx-\omega t)},\end{array}$$(A.21)

where f₊(v_x) and f₋(v_x) are not equal, according to what we show above. Introducing this fluctuation in the Vlasov equation we now obtain the following:

$$\begin{array}{cc}\hfill -i\omega e^{-i\omega t}(f_{+}{e}^{\mathit{ikx}}+f_{-}{e}^{-ikx})& +v_{x}ik{e}^{-i\omega t}(f_{+}{e}^{\mathit{ikx}}-f_{-}{e}^{-ikx})\hfill \\ \hfill & +\frac{q_s}{m_s}\delta G(x,t)\frac{\partial f_{0s}(v_x)}{\partial v_x}=0.\hfill \end{array}$$(A.22)

To have a proper solution, δG(x, t) must be of the form

$$\begin{array}{c}\hfill \delta G(x,t)=e^{-i\omega t}(G_{+}{e}^{\mathit{ikx}}+G_{-}{e}^{-ikx});\end{array}$$(A.23)

in this case the temporal phases cancel, and to have a solution the spatial phases are balanced in Eq. (A.22) if the following conditions are satisfied:

$$\begin{array}{cc}\hfill f_{+}& =-i\frac{q_s}{m_s}{G}_{+}\frac{1}{\omega -k{v}_x}\frac{\partial f_{0s}(v_x)}{\partial v_x},\hfill \\ \hfill f_{-}& =-i\frac{q_s}{m_s}{G}_{-}\frac{1}{\omega +k{v}_x}\frac{\partial f_{0s}(v_x)}{\partial v_x}.\hfill \end{array}$$(A.24)

These expressions provide the fluctuation associated with each travelling wave and are in agreement with Eqs. (A.6) and (A.17) when we performed the calculations for each propagating wave independently. Hereafter we assume that k is always positive. Now we can calculate the density fluctuation associated with the total distribution function based on Eq. (A.21) by performing the integration in velocities. The result is

$$\begin{array}{cc}\hfill \delta n_s(x,t)& =i{n}_{0s}\frac{q_s}{m_s}\frac{G_{+}}{k{w}_{0s}^2}{Z}^{\prime }({\xi }_s)e^{-i\omega t}{e}^{\mathit{ikx}}\hfill \\ \hfill & -i{n}_{0s}\frac{q_s}{m_s}\frac{G_{-}}{k{w}_{0s}^2}{Z}^{\prime }({\xi }_s)e^{-i\omega t}{e}^{-ikx}.\hfill \end{array}$$(A.25)

Interestingly, this expression simplifies the problem if G₊ = −G₋ since we obtain the simple form

$$\begin{array}{c}\hfill \delta n_s(x,t)=i{n}_{0s}\frac{q_s}{m_s}\frac{G_{+}}{k{w}_{0s}^2}{Z}^{\prime }({\xi }_s)e^{-i\omega t}(e^{\mathit{ikx}}+e^{-ikx}),\end{array}$$(A.26)

where we clearly have the superposition of the right- and left-propagating wave with the same amplitude. If we choose that G₊ = G₋, we simply have

$$\begin{array}{c}\hfill \delta n_s(x,t)=i{n}_{0s}\frac{q_s}{m_s}\frac{G_{+}}{k{w}_{0s}^2}{Z}^{\prime }({\xi }_s)e^{-i\omega t}(e^{\mathit{ikx}}-e^{-ikx}).\end{array}$$(A.27)

We note that the addition of exponentials in Eq. (A.26) reduces to a cosinus while the subtraction in Eq. (A.27) leads to a sinus. This is precisely the profile we obtained previously; see Eq. (16).

Appendix B: Forward- and backward-propagating waves

It is straightforward to get from Eq. (A.4) that

$$\begin{array}{c}\hfill \frac{\partial f_{0s}(v_x)}{\partial v_x}=\frac{n_{0s}}{\sqrt{\pi }{w}_{0s}}\frac{-2(v_x-u_{x0s})}{w_{0s}^2}{e}^{-\frac{{(v_x-u_{x0s})}^2}{w_{0s}^2}};\end{array}$$(B.1)

therefore, we have from Eq. (A.8) that for a wave propagating to the right

$$\begin{array}{c}\hfill \delta n_s(x,t)=i\frac{q_s}{m_s}\delta G_s(x,t)\frac{2n_{0s}}{\sqrt{\pi }{w}_{0s}^3}{I}_s(\omega ,k),\end{array}$$(B.2)

where

$$\begin{array}{c}\hfill I_s(\omega ,k)={\displaystyle {\int }_{-\infty }^{+\infty }\frac{(v_x-u_{x0s})e^{-\frac{{(v_x-u_{x0s})}^2}{w_{0s}^2}}}{\omega -k{v}_x}d{v}_x;}\end{array}$$(B.3)

by setting $\lambda =\frac{{\mathit{v}}_x-u_{x0s}}{{\mathit{w}}_{0s}}$ the integral I_s(ω, k) becomes

$$\begin{array}{c}\hfill I_s(\omega ,k)=-\frac{w_{0s}}{k}{\displaystyle {\int }_{-\infty }^{+\infty }\frac{\lambda e^{-{\lambda }^2}}{\lambda -{\xi }_s}d\lambda ,}\end{array}$$(B.4)

where ${\xi }_s=\frac{\omega -k{u}_{x0s}}{k{\mathit{w}}_{0s}}$, which reduces to

$$\begin{array}{c}\hfill I_s(\omega ,k)=\frac{\sqrt{\pi }{w}_{0s}}{2k}{Z}^{\prime }({\xi }_s),\end{array}$$(B.5)

where we define Z′(ξ_s) = − 2(1 + ξ_s Z(ξ_s)) expressed in terms of the plasma dispersion function Z(ξ_s) (e.g., the Fried and Conte function) defined as

$$\begin{array}{c}\hfill Z({\xi }_s)=\frac{1}{\sqrt{\pi }}{\displaystyle {\int }_{-\infty }^{+\infty }\frac{e^{-{\lambda }^2}}{\lambda -{\xi }_s}d\lambda .}\end{array}$$(B.6)

Finally, substituting Eq. (B.5) into Eq. (B.2) we obtain Eq. (A.9). The integral form of the Zeta function is in principle only applicable when Im(ξ_s) > 0.

For the left-propagating wave with wavenumber −k we have

$$\begin{array}{c}\hfill I_s(\omega ,-k)={\displaystyle {\int }_{-\infty }^{+\infty }\frac{v_x{e}^{-\frac{{(v_x)}^2}{w_{0s}^2}}}{\omega +k{v}_x}d{v}_x;}\end{array}$$(B.7)

by setting $\lambda =-\frac{{\mathit{v}}_x}{{\mathit{w}}_{0s}}$ the integral now becomes

$$\begin{array}{cc}\hfill I_s(\omega ,-k)& =\frac{w_{0s}}{k}{\displaystyle {\int }_{+\infty }^{-\infty }\frac{\lambda e^{-{\lambda }^2}}{{\xi }_s-\lambda }d\lambda =-\frac{w_{0s}}{k}{\int }_{-\infty }^{+\infty }\frac{\lambda e^{-{\lambda }^2}}{{\xi }_s-\lambda }d\lambda =}\hfill \\ \hfill & \frac{w_{0s}}{k}{\displaystyle {\int }_{-\infty }^{+\infty }\frac{\lambda e^{-{\lambda }^2}}{\lambda -{\xi }_s}d\lambda ,}\hfill \end{array}$$(B.8)

where ${\xi }_s=\frac{\omega }{k{\mathit{w}}_{0s}}$, and the integral reduces to

$$\begin{array}{c}\hfill I_s(\omega ,-k)=-\frac{\sqrt{\pi }{w}_{0s}}{2k}{Z}^{\prime }({\xi }_s).\end{array}$$(B.9)

This expression is valid if Im(ξ_s) > 0. We note the difference in the sign between Eq. (B.9) (left-propagating wave) and Eq. (B.5) (right-propagating wave). The previous analysis takes into account the effect of the sign of k in the dispersion relation, in agreement with the definition given in Podesta & Gary (2011) (see their Eq. (26)).

Appendix C: The standing Alfvén wave solution

Let us construct a standing Alfvénic wave solution using the results for circularly-propagating waves. The idea is to superimpose two waves with the same frequency and same wavenumber but travelling in opposite directions. This is easily done in MHD because of the degeneracy of Alfvén waves. In the multifluid approach used in this work changing ω by −ω leads to another solution of the dispersion relation. Nevertheless, the symmetry regarding the wavenumber is still present, as we can see from Eq. (48). Therefore, the propagation direction of the wave is easily reversed by changing the sign of k₀. We consider first right-handed Alfvén waves. The perturbed velocity and magnetic field components are of the form

$$\begin{array}{cc}\hfill u_z& =u_{0}cos(k_{0}x-{\omega }_{0}t),\hfill \end{array}$$(C.1)

$$\begin{array}{cc}\hfill B_z& =-u_0\frac{B_0}{v_{\mathrm{A}}^2}\frac{{\omega }_0}{k_0}cos(k_{0}x-{\omega }_{0}t),\hfill \end{array}$$(C.2)

$$\begin{array}{cc}\hfill u_y& =u_{0}sin(k_{0}x-{\omega }_{0}t),\hfill \end{array}$$(C.3)

$$\begin{array}{cc}\hfill B_y& =-u_0\frac{B_0}{v_{\mathrm{A}}^2}\frac{{\omega }_0}{k_0}sin(k_{0}x-{\omega }_{0}t),\hfill \end{array}$$(C.4)

where we use Eq. (50) to use the proper amplitude and sign of the magnetic field in terms of the initial velocity amplitude u₀. It is important to note that the amplitude of the perturbed magnetic field, $B_{\perp }=\sqrt{B_x^2+B_{\mathit{y}}^2}$ is a constant, independent of space and time (also true for the perpendicular velocity).

The corresponding backward wave (we change k₀ to −k₀) has the same frequency (see Eq. (48)) and using again Eq. (50) we have

$$\begin{array}{cc}\hfill u_z& =u_{0}cos(-k_{0}x-{\omega }_{0}t),\hfill \end{array}$$(C.5)

$$\begin{array}{cc}\hfill B_z& =u_0\frac{B_0}{v_{\mathrm{A}}^2}\frac{{\omega }_0}{k_0}cos(-k_{0}x-{\omega }_{0}t),\hfill \end{array}$$(C.6)

$$\begin{array}{cc}\hfill u_y& =u_{0}sin(-k_{0}x-{\omega }_{0}t),\hfill \end{array}$$(C.7)

$$\begin{array}{cc}\hfill B_y& =u_0\frac{B_0}{v_{\mathrm{A}}^2}\frac{{\omega }_0}{k_0}sin(-k_{0}x-{\omega }_{0}t).\hfill \end{array}$$(C.8)

Now we just need to make a linear combination of the two waves to obtain the standing solution. We choose to subtract the forward wave from the backward wave and using standard trigonometric formulas the result is written as

$$\begin{array}{cc}\hfill u_z& =u_{0}sin({\omega }_{0}t)sin(k_{0}x),\hfill \end{array}$$(C.9)

$$\begin{array}{cc}\hfill B_z& =-u_0\frac{B_0}{v_{\mathrm{A}}^2}\frac{{\omega }_0}{k_0}cos({\omega }_{0}t)cos(k_{0}x),\hfill \end{array}$$(C.10)

$$\begin{array}{cc}\hfill u_y& =u_{0}cos({\omega }_{0}t)sin(k_{0}x),\hfill \end{array}$$(C.11)

$$\begin{array}{cc}\hfill B_y& =u_0\frac{B_0}{v_{\mathrm{A}}^2}\frac{{\omega }_0}{k_0}sin({\omega }_{0}t)cos(k_{0}x),\hfill \end{array}$$(C.12)

where a common multiplicative constant has been eliminated. It is evident from the previous expressions that, contrary to the propagating case, B_⊥ is not constant now. To properly excite a standing right-hand wave in the simulations at t = 0 we just need to impose a perturbation of the following form:

$$\begin{array}{cc}\hfill u_z& =0,\hfill \end{array}$$(C.13)

$$\begin{array}{cc}\hfill B_z& =-u_0\frac{B_0}{v_{\mathrm{A}}^2}\frac{{\omega }_0}{k_0}cos(k_{0}x),\hfill \end{array}$$(C.14)

$$\begin{array}{cc}\hfill u_y& =u_{0}sin(k_{0}x),\hfill \end{array}$$(C.15)

$$\begin{array}{cc}\hfill B_y& =0.\hfill \end{array}$$(C.16)

A similar approach is used to derive the expression for the standing left-handed polarised wave. Now the dispersion relation is given by Eq. (48), and Eq. (50) does not change. For the forward wave we have

$$\begin{array}{cc}\hfill u_z& =-u_{0}cos(k_{0}x-{\omega }_{0}t),\hfill \end{array}$$(C.17)

$$\begin{array}{cc}\hfill B_z& =u_0\frac{B_0}{v_{\mathrm{A}}^2}\frac{{\omega }_0}{k_0}cos(k_{0}x-{\omega }_{0}t),\hfill \end{array}$$(C.18)

$$\begin{array}{cc}\hfill u_y& =u_{0}sin(k_{0}x-{\omega }_{0}t),\hfill \end{array}$$(C.19)

$$\begin{array}{cc}\hfill B_y& =-u_0\frac{B_0}{v_{\mathrm{A}}^2}\frac{{\omega }_0}{k_0}sin(k_{0}x-{\omega }_{0}t),\hfill \end{array}$$(C.20)

while for the backward wave (k₀ → −k₀) we have

$$\begin{array}{cc}\hfill u_z& =-u_{0}cos(-k_{0}x-{\omega }_{0}t),\hfill \end{array}$$(C.21)

$$\begin{array}{cc}\hfill B_z& =-u_0\frac{B_0}{v_{\mathrm{A}}^2}\frac{{\omega }_0}{k_0}cos(-k_{0}x-{\omega }_{0}t),\hfill \end{array}$$(C.22)

$$\begin{array}{cc}\hfill u_y& =u_{0}sin(-k_{0}x-{\omega }_{0}t),\hfill \end{array}$$(C.23)

$$\begin{array}{cc}\hfill B_y& =u_0\frac{B_0}{v_{\mathrm{A}}^2}\frac{{\omega }_0}{k_0}sin(-k_{0}x-{\omega }_{0}t).\hfill \end{array}$$(C.24)

Combining the previous waves we obtain

$$\begin{array}{cc}\hfill u_z& =u_{0}sin({\omega }_{0}t)sin(k_{0}x),\hfill \end{array}$$(C.25)

$$\begin{array}{cc}\hfill B_z& =-u_0\frac{B_0}{v_{\mathrm{A}}^2}\frac{{\omega }_0}{k_0}cos({\omega }_{0}t)cos(k_{0}x),\hfill \end{array}$$(C.26)

$$\begin{array}{cc}\hfill u_y& =-u_{0}cos({\omega }_{0}t)sin(k_{0}x),\hfill \end{array}$$(C.27)

$$\begin{array}{cc}\hfill B_y& =-u_0\frac{B_0}{v_{\mathrm{A}}^2}\frac{{\omega }_0}{k_0}sin({\omega }_{0}t)cos(k_{0}x).\hfill \end{array}$$(C.28)

Hence, to excite a circularly left-handed Alfvén wave at t = 0 we have to impose the following initial perturbation:

$$\begin{array}{cc}& u_z=0,\hfill \end{array}$$(C.29)

$$\begin{array}{cc}& B_z=-u_0\frac{B_0}{v_{\mathrm{A}}^2}\frac{{\omega }_0}{k_0}cos(k_{0}x),\hfill \end{array}$$(C.30)

$$\begin{array}{cc}& u_y=-u_{0}sin(k_{0}x),\hfill \end{array}$$(C.31)

$$\begin{array}{cc}& B_y=0.\hfill \end{array}$$(C.32)

Appendix D: Kinetic effects on ion-acoustic waves driven by Alfvén waves

We carry out a new regular perturbation method on the Vlasov equation to introduce kinetic effects by implementing the method described in Section 3.1. To begin the new analysis scheme we first introduce the perturbation expansion of the macroscopic fluid quantities:

$$\begin{array}{cc}\hfill \mathbf{u}& =ϵ{\mathbf{u}}_{\perp }+{ϵ}^2\delta u_x{\widehat{\mathbf{e}}}_x+{ϵ}^3{\mathbf{u}}_{\perp }^{\prime },\hfill \\ \hfill \mathbf{B}& =B_0{\widehat{\mathbf{e}}}_x+ϵ{\mathbf{B}}_{\perp }+{ϵ}^3{\mathbf{B}}_{\perp }^{\prime },\hfill \\ \hfill \mathbf{E}& =ϵ{\mathbf{E}}_{\perp }+{ϵ}^2\delta E_x{\widehat{\mathbf{e}}}_x+{ϵ}^3{\mathbf{E}}_{\perp }^{\prime },\hfill \\ \hfill n& =n_0+{ϵ}^2\delta n,\hfill \\ \hfill p& =p_0+{ϵ}^2\delta p.\hfill \end{array}$$(D.1)

The analysis begins with the linearisation scheme of the velocity distribution function

$$\begin{array}{cc}\hfill F_s(x,\mathbf{v},t)& =F_{0s}(\mathbf{v})+{ϵ}^2\delta F_s(x,\mathbf{v},t).\hfill \end{array}$$(D.2)

We note that the perturbation expansion of the velocity distribution does not contain first-order terms in ϵ. This is because in this order an incompressible transverse Alfvén wave is imposed, which has no density perturbation associated with it. Furthermore, the zero-order expansion is chosen to be independent of time and space in the following form:

$$\begin{array}{cc}\hfill F_{0s}(\mathbf{v})& =\frac{n_{0s}}{\sqrt{\pi }{w}_{0s}}\delta ({\mathbf{v}}_{\perp }-{{\mathbf{u}}_{\perp }}_{0s})e^{-\frac{{(v_x-u_{x0s})}^2}{w_{0s}^2}}.\hfill \end{array}$$(D.3)

In these equations the velocity vector v again corresponds to the particle distribution velocities, whereas u_⊥0s corresponds to the transverse Alfvén wave fluid velocity perturbation of the initial wave (given in Appendix C).

Inserting the perturbation expansion in Eqs. (D.1) and (D.2) into the Vlasov equation, we obtain a hierarchy of equations at different orders in the perturbation expansion of the distribution function. At the zeroth order in ϵ the Vlasov equation is satisfied exactly for the equilibrium state. To the first order in ϵ the Vlasov equation is also satisfied since there is no first-order velocity distribution in the expansion nor first-order density perturbation. At the second order in ϵ the Vlasov equation becomes

$$\begin{array}{cc}\hfill \frac{\partial \delta F_s(x,\mathbf{v},t)}{\partial t}& +\mathbf{v}·\mathrm{\nabla }\delta F_s(x,\mathbf{v},t)\hfill \\ \hfill & +\frac{q_s}{m_s}{(\mathbf{E}(x,t)+\frac{\mathbf{v}\times \mathbf{B}(x,t)}{c})}_x^{(0)}·\frac{\partial \delta F_s(x,\mathbf{v},t)}{\partial {\mathbf{v}}_x}\hfill \\ \hfill & +\frac{q_s}{m_s}{(\mathbf{E}(x,t)+\frac{\mathbf{v}\times \mathbf{B}(x,t)}{c})}_{\perp }^{(0)}·\frac{\partial \delta F_s(x,\mathbf{v},t)}{\partial {\mathbf{v}}_{\perp }}\hfill \\ \hfill & +\frac{q_s}{m_s}{(\mathbf{E}(x,t)+\frac{\mathbf{v}\times \mathbf{B}(x,t)}{c})}_x^{(2)}·\frac{\partial F_{0s}(x,\mathbf{v},t)}{\partial {\mathbf{v}}_x}\hfill \\ \hfill & +\frac{q_s}{m_s}{(\mathbf{E}(x,t)+\frac{\mathbf{v}\times \mathbf{B}(x,t)}{c})}_{\perp }^{(2)}·\frac{\partial F_{0s}(x,\mathbf{v},t)}{\partial {\mathbf{v}}_{\perp }}=0,\hfill \end{array}$$(D.4)

where the force terms have been separated into parallel (x) and perpendicular ( ⊥ ) components relative to the mean magnetic field direction and the superscripts in brackets, (0) and (2), represent the zeroth- and second-order force terms defined below. Thus, to the second order in ϵ the Vlasov equation reduces to

$$\begin{array}{cc}\hfill \frac{\partial \delta F_s(x,\mathbf{v},t)}{\partial t}& +\mathbf{v}·\mathrm{\nabla }\delta F_s(x,\mathbf{v},t)\hfill \\ \hfill =& -\frac{q_s}{m_s}(\delta E_x(x,t)+\frac{{({\mathbf{u}}_{\perp }(x,t)\times {\mathbf{B}}_{\perp }(x,t))}_x}{c})·\frac{\partial F_{0s}(\mathbf{v})}{\partial {\mathbf{v}}_x}·\hfill \end{array}$$(D.5)

At this point, in order to get a solution for the second-order distribution δF_s(x, v, t) we have to assume either a propagating or a standing solution. Equation (D.5) is equivalent to Eq. (A.5) indicating that it will eventually lead to the definition of the complex γ.

All Figures