ALFVÉN WAVES IN A PARTIALLY IONIZED TWO-FLUID PLASMA

Alfvén waves are a particular class of magnetohydrodynamic (MHD) waves driven by magnetic tension (Alfvén 1942). In a uniform and infinite plasma the motions of Alfvén waves are incompressible and polarized perpendicularly to the direction of the magnetic field (see, e.g., Hasegawa & Uberoi 1982; Cramer 2001; Goossens 2003). Alfvén waves are found in both laboratory and astrophysical plasmas (see review by Gekelman et al. 2011).

Since the pioneering works by, e.g., Piddington (1956) and Kulsrud & Pearce (1969), it is known that partial ionization of the plasma affects the dynamics of Alfvén waves. The feature most extensively investigated in the literature is the wave damping due to collisions between ions and neutrals, although other effects as, e.g., the existence of cutoff values of the wavelength are also an important consequence of partial ionization (Kulsrud & Pearce 1969).

Most of the works that studied Alfvén waves in partially ionized plasmas adopted the so-called single-fluid approximation (see, e.g., Braginskii 1965). The single-fluid approximation assumes a strong coupling between ions and neutrals. For an MHD wave, this condition means that the wave frequency has to be much lower than the frequency at which ions and neutrals collide. In other words, it is necessary that in one period there are enough collisions for ions and neutrals to behave as one fluid. This restriction is fulfilled in, e.g., the partially ionized solar plasma and so the single-fluid approximation is usually adopted in that case (see, e.g., De Pontieu et al. 2001; Khodachenko et al. 2004; Forteza et al. 2007; Soler et al. 2009, among others).

An alternative approach is the multi-fluid theory (see, e.g., Zaqarashvili et al. 2011b), which considers the various species in the plasma as separate fluids. In the multi-fluid description no restriction is imposed on the relative values of the wave frequency and the collision frequency, although the mathematical treatment is usually more complicated than in the single-fluid case. However, the multi-fluid theory has the advantage that it is more general than the single-fluid approximation and can be used regardless of the value of the collision frequency. Hence, the multi-fluid theory is adequate to study MHD waves in those situations where the single-fluid approximation does not apply. This may be the case of molecular clouds (see, e.g., Pudritz 1990; Balsara 1996; Mouschovias et al. 2011). A particular form of the multi-fluid theory is the two-fluid theory in which ions and electrons are considered together as an ion–electron fluid, while neutrals form another fluid that interacts with the ion–electron fluid by means of collisions. For the investigation of MHD waves this approach was followed by, e.g., Kumar & Roberts (2003), Zaqarashvili et al. (2011b), Mouschovias et al. (2011), and Soler et al. (2012).

Despite the existing literature on this topic (see the references in the above paragraphs), the purpose of the present article is to revisit the theoretical investigation of Alfvén waves in partially ionized plasmas using the two-fluid theory. Our reasons for tackling this task are the following.

First of all, the existing papers in the literature often focus on very specific situations as, e.g., the interstellar medium (e.g., Kulsrud & Pearce 1969), molecular clouds (see, e.g., Pudritz 1990; Balsara 1996; Mouschovias et al. 2011), and solar plasmas (Kumar & Roberts 2003; Zaqarashvili et al. 2011b; Soler et al. 2012) among other cases. Here our aim is to keep the investigation as general as possible. To do so we take the neutral–ion collision frequency and the plasma ionization degree as free parameters. This makes the results of the present article widely applicable. We put emphasis on the mathematical transparency and on the finding of approximate analytic solutions.

In addition, there has been recently some confusion about the existence of cutoff wavelengths for Alfvén waves in a partially ionized plasma (see Zaqarashvili et al. 2011b, 2012). While Zaqarashvili et al. (2012) have shown that the presence of cutoffs in the single-fluid approximation is a mathematical artifact, the existence of cutoffs in the two-fluid case is a real physical phenomenon (e.g., Kulsrud & Pearce 1969; Pudritz 1990; Kamaya & Nishi 1998; Mouschovias et al. 2011). An important goal of the present article is to stress the existence of physical cutoffs for Alfvén waves in a two-fluid plasma.

Finally, unlike previous works that are restricted to the normal mode analysis, here we combine results of normal modes with the solution of the initial-value problem. This procedure allows us to investigate how strong is the coupling between the perturbations in the ionized fluid and the neutral fluid depending on the relative values of the wave frequency and the neutral–ion collision frequency.

This paper is organized as follows. Section 2 contains the description of the equilibrium configuration and the basic equations of the two-fluid theory. Alfvén waves are investigated following a normal mode analysis in Section 3, while the initial-value problem is solved in Section 4. Later, Section 5 contains an application to the solar atmospheric plasma, and Section 6 discusses the implications of partial ionization for the energy flux of Alfvén waves. Finally, the conclusions of this work are given in Section 7.

We consider a partially ionized medium composed of ions, electrons, and neutrals. We use the two-fluid theory in which ions and electrons are considered together as an ion–electron fluid, i.e., the ionized fluid, while neutrals form another fluid that interacts with the ionized fluid by means of collisions (see, e.g., Zaqarashvili et al. 2011b; Soler et al. 2012). In all the following expressions, the subscripts "i" and "n" refer to the ionized fluid and the neutral fluid, respectively.

The equilibrium is made of a uniform and unbounded partially ionized plasma. We use Cartesian coordinates. The equilibrium magnetic field is straight and constant along the z-direction, namely ${\bf B}=B\, \hat{z}$. We ignore the effect of gravity. We also assume that the equilibrium is static so that there are no equilibrium flows. The governing equations for the various species composing the plasma can be found in, e.g., Zaqarashvili et al. (2011b). Here we restrict ourselves to the study of linear perturbations superimposed on the equilibrium state. Hence, the governing equations are linearized. The resulting equations are

$$\begin{eqnarray} \rho _{\rm i}\frac{\partial {\bf v}_{\rm i}}{\partial t} &=& {-} \nabla p_{\rm i} + \frac{1}{\mu }\left(\nabla \times {\bf b} \right) \times {\bf B} -\alpha _{\rm in}\left({\bf v}_{\rm i} - {\bf v}_{\rm n}\right), \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \rho _{\rm n}\frac{\partial {\bf v}_{\rm n}}{\partial t} &=& {-} \nabla p_{\rm n} -\alpha _{\rm in}\left({\bf v}_{\rm n} - {\bf v}_{\rm i}\right), \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \frac{\partial {\bf b}}{\partial t} &=& \nabla \times \left({\bf v}_{\rm i} \times {\bf B} \right), \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \frac{\partial p_{\rm i}}{\partial t} &=& {-} \gamma P_{\rm i} \nabla \cdot {\bf v}_{\rm i}, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} \frac{\partial p_{\rm n}}{\partial t} &=& {-} \gamma P_{\rm n} \nabla \cdot {\bf v}_{\rm n}, \end{eqnarray} \tag{ 5 }$$

where v_i, p_i, P_i, and ρ_i are the velocity perturbation, pressure perturbation, equilibrium pressure, and equilibrium density of the ionized fluid, v_n, p_n, P_n, and ρ_n are the respective quantities but for the neutral fluid, b is the magnetic field perturbation, μ is the magnetic permeability, γ is the adiabatic index, and α_in is the ion–neutral friction coefficient. In the specific case of a hydrogen plasma, the expression of α_in is given by Braginskii (1965), namely

$$\begin{equation} \alpha _{\rm in}= \frac{1}{2} \frac{\rho _{\rm i}\rho _{\rm n}}{m_{\rm n}}\sqrt{\frac{16 k_{\rm B} T}{\pi m_{\rm i}}}\sigma _{\rm in}, \end{equation} \tag{ 6 }$$

where m_i and m_n are the ion and neutral masses, respectively (m_i ≈ m_n for hydrogen), k_B is Boltzmann's constant, T is the plasma temperature, and σ_in is the collision cross section. In the following analysis we do not use this expression of α_in since we take α_in as a free parameter. We do so to conveniently control the strength of the ion–neutral friction force. Equation (6) is used in the application to solar plasmas done in Section 5.

We perform a Fourier analysis of the perturbations in space. In linear theory, an arbitrary perturbation can be represented by the superposition of Fourier components, so that we can restrict ourselves to study particular Fourier components. Therefore, the spatial dependence of perturbations is put proportional to exp (ik_xx + ik_yy + ik_zz), where k_x, k_y, and k_z are the components of the wavenumber in the x-, y-, and z-directions, respectively. In a uniform and infinite plasma, Alfvén waves are the only MHD modes that propagate vorticity perturbations (see, e.g., Cramer 2001; Goossens 2003). In addition, Alfvén waves are incompressible and their motions are confined to perpendicular planes to the magnetic field, i.e., v_{i, z} = v_{n, z} = 0. Therefore, an appropriate quantity to describe Alfvén waves is the vorticity component along the magnetic field direction. By working with vorticity perturbations we are able to decouple Alfvén waves from magnetoacoustic waves. We define Γ_i and Γ_n as the z-components of vorticity of the ionized fluid and the neutral fluid, respectively,

$$\begin{eqnarray} \Gamma _{\rm i}&=& \left(\nabla \times {\bf v}_{\rm i} \right) \cdot \hat{z} = i k_x v_{{\rm i}, y} - i k_y v_{{\rm i}, x}, \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} \Gamma _{\rm n}&=& \left(\nabla \times {\bf v}_{\rm n} \right) \cdot \hat{z} = i k_x v_{{\rm n},y} - i k_y v_{{\rm n},x}. \end{eqnarray} \tag{ 8 }$$

Note that in the reference frame in which k_y = 0, Γ_i and Γ_n are proportional to v_{i, y} and v_{n, y}, respectively. We combine Equations (1)–(5) and after some algebraic manipulations, we obtain the two following equations involving Γ_i and Γ_n only, namely

$$\begin{eqnarray} \rho _{\rm i}\frac{\partial ^2 \Gamma _{\rm i}}{\partial t^2} + \alpha _{\rm in}\frac{\partial \Gamma _{\rm i}}{\partial t} + \rho _{\rm i}k_z^2 c_{\mathrm{A}}^2 \Gamma _{\rm i}&=& \alpha _{\rm in}\frac{\partial \Gamma _{\rm n}}{\partial t}, \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} \rho _{\rm n}\frac{\partial \Gamma _{\rm n}}{\partial t} + \alpha _{\rm in}\Gamma _{\rm n}&=& \alpha _{\rm in}\Gamma _{\rm i}, \end{eqnarray} \tag{ 10 }$$

where $c_{\mathrm{A}}= B/\sqrt{\mu \rho _{\rm i}}$ is the Alfvén velocity. Note that the Alfvén velocity is here defined using the density of the ionized fluid only. Equations (9) and (10) are the governing equations for linear vorticity perturbations and, therefore, they are the governing equations of Alfvén waves.

For the subsequent analysis we define the ionization fraction, χ, the ion–neutral collision frequency, ν_in, and the neutral–ion collision frequency, ν_ni, as follows:

$$\begin{equation} \chi = \frac{\rho _{\rm n}}{\rho _{\rm i}}, \quad \nu _{\rm in} = \frac{\alpha _{\rm in}}{\rho _{\rm i}}, \quad \nu _{\rm ni} = \frac{\alpha _{\rm in}}{\rho _{\rm n}}. \end{equation} \tag{ 11 }$$

Since the collision frequencies are related by ρ_iν_in = ρ_nν_ni, we use ν_ni in all the following expressions for simplicity. Note that when ρ_i ≠ ρ_n, ν_in ≠ ν_ni, meaning that the ion–neutral and neutral–ion collision frequencies are different (see a discussion on this issue in Zaqarashvili et al. 2011a).

Here we perform a normal mode analysis. The temporal dependence of the perturbations is put proportional to exp (− iωt), where ω is the angular frequency. From Equation (10), we express Γ_n in terms of Γ_i and insert the expression in Equation (9). We arrive at an equation involving Γ_i only, namely

$$\begin{equation} \mathcal {D}\left(\omega, k_z \right) \Gamma _{\rm i}= 0, \end{equation} \tag{ 12 }$$

with

$$\begin{equation} \mathcal {D}\left(\omega, k_z \right) = \omega ^3 + i \left(1+\chi \right) \nu _{\mathrm{ni}}\omega ^2 - k_z^2 c_{\rm A}^2 \omega - i \nu _{\mathrm{ni}}k_z^2 c_{\rm A}^2. \end{equation} \tag{ 13 }$$

For Γ_i ≠ 0, the solutions to Equation (12) must satisfy $\mathcal {D}\left(\omega, k_z \right) = 0$, i.e.,

$$\begin{equation} \omega ^3 + i \left(1+\chi \right) \nu _{\mathrm{ni}}\omega ^2 - k_z^2 c_{\rm A}^2 \omega - i \nu _{\mathrm{ni}}k_z^2 c_{\rm A}^2=0. \end{equation} \tag{ 14 }$$

Equation (14) is the dispersion relation of Alfvén waves. Although with different notations, Equation (14) is equivalent to the dispersion relations previously found by, e.g., Piddington (1956), Kulsrud & Pearce (1969), Pudritz (1990), Martin et al. (1997), Kamaya & Nishi (1998), Kumar & Roberts (2003), Zaqarashvili et al. (2011b), and Mouschovias et al. (2011). In the absence of collisions, ν_ni = 0 and Equation (14) becomes

$$\begin{equation} \omega \left(\omega ^2-k_z^2c_{\mathrm{A}}^2 \right) = 0. \end{equation} \tag{ 15 }$$

From Equation (15) we get the classic dispersion relation of Alfvén waves in an ideal plasma, namely $\omega ^2=k_z^2c_{\mathrm{A}}^2$, and an additional mode with ω = 0. The general situation ν_ni ≠ 0 is investigated next.

3.1. Standing Waves

We focus first on standing waves. Hence we assume a real wavenumber, k_z, and solve the dispersion relation (Equation (14)) to obtain the complex frequency, ω = ω_R + iω_I, with ω_R and ω_I the real and imaginary parts of ω, respectively. Since ω is complex the amplitude of perturbations is multiplied by the factor exp (ω_It), with ω_I < 0. Therefore, the perturbations are damped in time.

Equation (14) is a cubic equation so it has three solutions. Unfortunately, the exact analytic solution to Equation (14) is too complicated to shed any light on the physics. However, we can investigate the nature of the solutions using the concept of the polynomial discriminant. We perform the change of variable ω = −is, so that Equation (14) becomes

$$\begin{equation} s^3 + \left(1+\chi \right) \nu _{\mathrm{ni}}s^2 + k_z^2 c_{\rm A}^2 s + \nu _{\mathrm{ni}}k_z^2 c_{\rm A}^2 = 0. \end{equation} \tag{ 16 }$$

Equation (16) is a cubic equation and all its coefficients are real. From Equation (16) we compute the discriminant, Λ, namely (see, e.g., Cohen 2000)

$$\begin{eqnarray} \Lambda = &-& k_z^2 c_{\rm A}^2 \left[ 4\left(1+\chi \right)^3 \nu _{\mathrm{ni}}^4 \right. \nonumber \\ &-& \left. (\chi ^2 + 20\chi -8)\nu _{\mathrm{ni}}^2 k_z^2 c_{\rm A}^2 +4 k_z^4 c_{\rm A}^4\right]. \end{eqnarray} \tag{ 17 }$$

The discriminant, Λ, is defined so that (1) Equation (16) has one real zero and two complex conjugate zeros when Λ < 0, (2) Equation (16) has a multiple zero and all the zeros are real when Λ = 0, and (3) Equation (16) has three distinct real zeros when Λ > 0. This classification is very relevant because the complex zeros of Equation (16) result in damped oscillatory solutions of Equation (14), whereas the real zeros of Equation (16) correspond to evanescent solutions of Equation (14).

It is instructive to consider again the paradigmatic situation in which there are no collisions between the two fluids, so we set ν_ni = 0. The discriminant becomes $\Lambda = -4k_z^6 c_{\rm A}^6<0$, which means that Equation (16) has one real zero and two complex conjugate zeros. Indeed, when ν_ni = 0 the zeros of Equation (16) are

$$\begin{equation} s = \pm i k_z c_{\mathrm{A}}, \quad s=0, \end{equation} \tag{ 18 }$$

which correspond to the following values of ω:

$$\begin{equation} \omega = \pm k_zc_{\mathrm{A}}, \qquad \omega =0. \end{equation} \tag{ 19 }$$

The two non-zero solutions correspond to the ideal Alfvén frequency, as expected.

We go back to the general case ν_ni ≠ 0. To determine the location where the nature of the solutions changes we set Λ = 0 and find the corresponding relation between the various parameters. For given ν_ni and χ, we find two different values of k_z, denoted by $k_z^+$ and $k_z^-$, which satisfy Λ = 0, namely

$$\begin{equation} k_z^\pm = \frac{\nu _{\rm ni}}{c_{\rm A}}\left[\frac{\chi ^2+20\chi -8}{8\left(1+\chi \right)^3} \pm \frac{\chi ^{1/2} \left(\chi -8 \right)^{3/2}}{8\left(1+\chi \right)^3}\right]^{-1/2}. \end{equation} \tag{ 20 }$$

Since k_z must be real, Equation (20) imposes a condition on the minimum value of χ which allows Λ = 0. This minimum value is χ = 8 and the corresponding critical k_z is $k_z^+ = k_z^- = 3\sqrt{3} \nu _{\mathrm{ni}}/c_{\mathrm{A}}$. When χ > 8, Equation (20) gives $k_z^+ < k_z^-$. For k_z outside the interval $(k_z^+,k_z^-),$ we have Λ < 0 so that there are two propagating Alfvén waves and one evanescent solution. For $k_z\in (k_z^+,k_z^-)$ we have Λ > 0 and so all three zeros of Equation (16) are real, i.e., they correspond to purely imaginary solutions of Equation (14). There is no propagation of Alfvén waves for $k_z\in (k_z^+,k_z^-)$. We call this interval the cutoff region. To the best of our knowledge, Kulsrud & Pearce (1969) were the first to report on the existence of a cutoff region of wavenumbers for Alfvén waves in a partially ionized two-fluid plasma, while studying the propagation of cosmic rays. These cutoffs also appear in the works by, e.g., Pudritz (1990), Kumar & Roberts (2003), and Mouschovias et al. (2011). The cutoff wavenumbers were ignored by Zaqarashvili et al. (2011b), who stated that there is always a solution of Equation (14) with a non-zero real part. Zaqarashvili et al. (2011b) probably reached this wrong conclusion because they never took χ > 8 in their computations. Here we clearly see that the three solutions of Equation (14) are purely imaginary when χ > 8 and $k_z\in (k_z^+,k_z^-)$.

Kamaya & Nishi (1998) discussed the physical reason for the existence of a range of cutoff wavenumbers in weakly ionized plasmas (see also Mouschovias 1987). When $k_z > k_z^-$ magnetic tension drives ions to oscillate almost freely, since the friction force is not strong enough to transfer significant inertia to neutrals. In this case, disturbances in the magnetic field affect only the ionized fluid as happens for classic Alfvén waves in fully ionized plasmas. Conversely, when $k_z < k_z^+$ the ion–neutral friction is efficient enough for neutrals to be nearly frozen into the magnetic field. After a perturbation, neutrals are dragged by ions almost instantly and both species oscillate together as a single fluid. The intermediate situation occurs when $k_z\in (k_z^+,k_z^-)$. In this case, a disturbance in the magnetic field decays due to friction before the ion–neutral coupling has had time to transfer the restoring properties of magnetic tension to the neutral fluid. In other words, neutral–ion collisions are efficient enough to dissipate perturbations in the magnetic field but, on the contrary, they are not efficient enough to transfer significant inertia to neutrals before the magnetic field perturbations have decayed. Hence, oscillations of the magnetic field are suppressed when $k_z\in (k_z^+,k_z^-)$. Additional insight on the physical behavior of the perturbations near the cutoff region is given in Section 3.1.2 by analyzing the forces acting on the fluids.

To avoid confusion we must inform the reader that the existence of cutoff wavenumbers of Alfvén waves discussed above is a purely two-fluid effect. These cutoff wavenumbers are not the cutoffs obtained in the single-fluid approximation (see, e.g., Balsara 1996; Forteza et al. 2007; Soler et al. 2009; Barceló et al. 2011). Zaqarashvili et al. (2012) have shown that the cutoff wavenumbers found in the single-fluid approximation are a mathematical artifact, i.e., they are caused by the approximations made when proceeding from the multi-fluid equations to single-fluid equations and are not connected to any real physical process. On the contrary, the cutoff wavenumbers found in the two-fluid case are physically and mathematically real and are caused by the two-fluid interaction between ions and neutrals (see Pudritz 1990; Kamaya & Nishi 1998; Mouschovias et al. 2011).

3.1.1. Approximate Analytic Solutions

We look for approximate analytic solutions to Equation (14) corresponding to standing modes. We assume that k_z is outside the cutoff interval $(k_z^+,k_z^-)$, so that Equation (14) has two complex solutions and one purely imaginary solution. This is the most interesting situation for the study of standing Alfvén waves since no oscillatory modes exist when k_z is within the cutoff interval. First, we look for an approximate expression for the two oscillatory solutions. To do so we write ω = ω_R + iω_I and insert this expression in Equation (14). We assume |ω_I| ≪ |ω_R| and neglect terms with $\omega _{\rm I}^2$ and higher powers. Hence, it is crucial for the validity of this approximation that k_z is not within or close to the cutoff region where ω_R = 0. After some algebraic manipulations we derive approximate expressions for ω_R and ω_I. For simplicity we omit the intermediate steps and give the final expressions, namely

$$\begin{eqnarray} \omega _{\rm R} &\approx & {\pm }k_z c_{\mathrm{A}}\sqrt{\frac{k_z^2 c_{\rm A}^2+\left(1 +\chi \right)\nu _{\mathrm{ni}}^2}{ k_z^2 c_{\rm A}^2+\left(1 +\chi \right)^2\nu _{\mathrm{ni}}^2}}, \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} \omega _{\rm I}&\approx & {-} \frac{\chi \nu _{\mathrm{ni}}}{2\left[ k_z^2 c_{\rm A}^2 +\left(1+\chi \right)^2 \nu _{\mathrm{ni}}^2\right]} k_z^2 c_{\rm A}^2. \end{eqnarray} \tag{ 22 }$$

On the other hand, the remaining purely imaginary, i.e., evanescent, solution is ω = i, with the approximation to given by

$$\begin{equation} \epsilon \approx - \nu _{\mathrm{ni}}\frac{k_z^2 c_{\rm A}^2+\left(1 +\chi \right)^2\nu _{\mathrm{ni}}^2}{k_z^2 c_{\rm A}^2+\left(1 +\chi \right)\nu _{\mathrm{ni}}^2}. \end{equation} \tag{ 23 }$$

When ν_ni = 0, we find ω_R = ±k_zc_A, ω_I = 0, and = 0; hence, we recover the solutions in the uncoupled case (Equation (19)).

It is useful to investigate the behavior of the solutions in the various limits of ν_ni. First, we consider the limit ν_ni ≪ k_zc_A, i.e., the case of low collision frequency, which means that the coupling between fluids is weak. Equations (21) and (23) simplify to

$$\begin{eqnarray} \omega _{\rm R} &\approx & {\pm } k_z c_{\mathrm{A}}, \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} \omega _{\rm I}&\approx & {-}\frac{\chi \nu _{\mathrm{ni}}}{2}, \end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray} \epsilon & \approx & {-} \nu _{\mathrm{ni}}. \end{eqnarray} \tag{ 26 }$$

In this limit ω_R coincides with its value in the ideal, uncoupled case and ω_I is independent of k_z. Hence, the damping of Alfvén waves does not depend on the wavenumber.

On the other hand, when ν_ni ≫ k_zc_A, i.e., the case of strong coupling between fluids, we find

$$\begin{eqnarray} \omega _{\rm R} &\approx & {\pm }\frac{k_z c_{\mathrm{A}}}{\sqrt{1+\chi }}, \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} \omega _{\rm I}&\approx & {-} \frac{\chi }{2\left(1+\chi \right)^2 } \frac{k_z^2 c_{\rm A}^2 }{\nu _{\mathrm{ni}}}, \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} \epsilon & \approx & {-} (1+\chi)\nu _{\mathrm{ni}}. \end{eqnarray} \tag{ 29 }$$

Now the expression of ω_R involves the factor $\sqrt{1+\chi }$ in the denominator, so that the larger the amount of neutrals, the lower ω_R compared to the value in the fully ionized case (see also Kumar & Roberts 2003; Soler et al. 2012). Now ω_I is proportional to $k_z^2$, meaning that the shorter the wavelength, the more efficient the damping.

3.1.2. Comparison with Numerical Results

Here we solve the full dispersion relation (Equation (14)) numerically and compare the numerical solutions with the previous approximations (Equations (21)–(23)). First, we set χ = 2 and vary the ratio ν_ni/k_zc_A between 10⁻² and 10². We compute ω_R/k_zc_A and ω_I/k_zc_A (see Figure 1). The agreement between numerical and analytic results is very good. There is no cutoff region for this choice of parameters because we have taken χ < 8. Regarding the real part of the frequency, we obtain that the oscillatory modes have ω_R/k_zc_A ≈ ±1 when ν_in/k_zc_A ≪ 1. When the ratio ν_ni/k_zc_A increases, ω_R/k_zc_A decreases until the value $\omega _{\rm R}/ k_z c_{\mathrm{A}}\approx \pm 1/\sqrt{1+\chi }$ is reached. This behavior is consistent with the analytic Equation (21). The evanescent mode has ω_R = 0 regardless the value of ν_ni/k_zc_A. The imaginary part of the frequency of the oscillatory modes tends to zero in both limits ν_ni/k_zc_A ≪ 1 and ν_ni/k_zc_A ≫ 1, while the damping is most efficient when ν_ni/k_zc_A ∼ 1. This result is also consistent with the analytic Equation (22), although the approximation underestimates the actual damping rate when ν_ni/k_zc_A ∼ 1. This discrepancy is a consequence of the weak damping approximation, which assumes |ω_I| ≪ |ω_R|. However, when ν_ni/k_zc_A ∼ 1, the numerical results show that |ω_I| and |ω_R| are of the same order; hence, the damping is strong. The imaginary part of the frequency of the evanescent mode is very well approximated by Equation (23).

Zoom In Zoom Out Reset image size

Next we increase ionization ratio to χ = 20 and compute the same results as before (Figure 2). Now there is a cutoff region because χ > 8. The cutoff region is correctly described by Equation (20). At the cutoff the oscillatory and evanescent modes interact and they become three purely imaginary solutions. Propagation is forbidden in this interval. As expected, the agreement between numerical and analytic results is not good near the cutoff region, but both results are in reasonably agreement far from the cutoff interval.

Zoom In Zoom Out Reset image size

To explore the physical behavior of the solutions near the cutoff region, we rewrite the momentum equations of ions (Equation (1)) and neutrals (Equation (2)) in the following forms:

$$\begin{eqnarray} \rho _{\rm i}\frac{\partial {\bf v}_{\rm i}}{\partial t} &=& {\bf T} - {\bf R}, \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} \rho _{\rm n}\frac{\partial {\bf v}_{\rm n}}{\partial t} &=& {\bf R}, \end{eqnarray} \tag{ 31 }$$

where T and R are the magnetic tension force and the friction force, respectively, given by

$$\begin{eqnarray} {\bf T} &=& {-} i \rho _{\rm i}\frac{k_z^2 c_{\mathrm{A}}^2}{\omega } {\bf v}_{\rm i}, \end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray} {\bf R} &=& \rho _{\rm n}\frac{\nu _{\mathrm{ni}}\omega }{\omega + i\nu _{\mathrm{ni}}} {\bf v}_{\rm i}. \end{eqnarray} \tag{ 33 }$$

In Equations (30) and (31) we have not included magnetic pressure and gas pressure forces because they do not affect Alfvén waves. Now we use the numerically obtained solutions for χ = 20 (Figure 2) to compute the moduli of T and R, namely ||T|| and ||R||, as functions of ν_ni/k_zc_A. Figure 3 displays the ratio ||T||/||R|| versus ν_ni/k_zc_A near the cutoff region. We have selected some locations in Figure 3, denoted by letters from a to e, to support the following discussion on the importance of the two forces.

Zoom In Zoom Out Reset image size

We start by analyzing the solutions on the left-hand side to the cutoff region. There are an oscillatory solution, a, and an evanescent solution, b. We find that ||T|| ≫ ||R|| for the oscillatory solution a, so that there is a net restoring force for ions in Equation (30). Magnetic tension is the dominant force and drives ions to oscillate almost freely, whereas neutrals are only slightly perturbed by the weak friction force in Equation (31). In the case of the evanescent solution b we obtain that ||T|| ≈ ||R||. This means that there is no net force acting on ions. The evanescent solution b only produces perturbations in the neutral fluid.

We turn to location c in Figure 3, i.e., within the cutoff interval. Here all the solutions are evanescent. The ratio ||T||/||R|| of the solution that was previously oscillatory decreases and becomes ||T||/||R|| < 1 before reaching the cutoff region. Now friction is the dominant force. Friction acts very efficiently in dissipating perturbations in the plasma before ions (and neutrals indirectly) have had time to feel the restoring force of magnetic tension. As a consequence, oscillatory modes are suppressed.

Finally, we analyze the forces acting on the solutions on the right-hand side to the cutoff region. Again, there are an oscillatory solution, d, and an evanescent solution, e. As happened for the oscillatory solution a, we find that ||T|| > ||R|| for the oscillatory solution d, although in this case the friction force is not negligible. Magnetic tension provides now the necessary restoring force for the oscillations of ions, while the friction force is responsible for dragging neutrals when ions move. Hence, both species tend to oscillate together. Solution d represents a collective oscillation of the whole plasma. On the contrary, magnetic tension is negligible for the evanescent solution, e. This mode is governed by the friction force alone and simply causes the decay of perturbations.

3.2. Propagating Waves

We move to the study of propagating waves. In the ideal fully ionized case the study of propagating waves is equivalent to that of standing waves. Here we shall see that ion–neutral collisions break this equivalence and propagating waves are worth being studied separately. For propagating waves we assume a real ω and solve the dispersion relation (Equation (14)) to find the complex wavenumber, k_z = k_z, R + ik_z, I, with k_{z, R} and k_{z, I} the real and imaginary parts of k_z, respectively. Then the amplitude of perturbations is multiplied by the factor exp (− k_{z, I}z), so that the perturbations are spatially damped. For real and positive ω, k_{z, R} > 0 corresponds to waves propagating toward the positive z-direction. Conversely, k_{z, R} < 0 corresponds to waves propagating toward the negative z-direction. In both situations, the sign of k_{z, I} is the same as that of k_{z, R}.

Equation (14) is a quadratic equation in k_z. The solution to Equation (14) is

$$\begin{equation} k_z^2 = \frac{\omega ^2}{c_{\mathrm{A}}^2} \frac{\omega + i (1+\chi)\nu _{\mathrm{ni}}}{\omega +i\nu _{\mathrm{ni}}}. \end{equation} \tag{ 34 }$$

When ν_ni = 0 we recover the ideal result $k_z^2 = \omega ^2/c_{\mathrm{A}}^2$. For ν_ni ≠ 0 we write k_z = k_z, R + ik_z, I and insert this expression in Equation (34). Then, it is possible to obtain the exact expression of $k_{z,\rm R}^2$, namely

$$\begin{eqnarray} k_{z,\rm R}^2 &=& \frac{1}{2} \frac{\omega ^2}{c_{\mathrm{A}}^2} \frac{\omega ^2 + (1+\chi)\nu _{\mathrm{ni}}^2}{\omega ^2 + \nu _{\mathrm{ni}}^2} \nonumber \\ &&\times \left[ 1 + \left(1 + \frac{\chi ^2\nu _{\mathrm{ni}}^2\omega ^2}{\left(\omega ^2 + (1+\chi)\nu _{\mathrm{ni}}^2\right)^2} \right)^{1/2} \right], \end{eqnarray} \tag{ 35 }$$

while the exact expression of $k_{z,\rm I}^2$ is

$$\begin{equation} k_{z,\rm I}^2 = k_{z,\rm R}^2 - \frac{\omega ^2}{c_{\mathrm{A}}^2} \frac{\omega ^2 + (1+\chi)\nu _{\mathrm{ni}}^2}{\omega ^2 + \nu _{\mathrm{ni}}^2}. \end{equation} \tag{ 36 }$$

Contrary to the case of standing waves there is no cutoff region for propagating Alfvén waves. It is always found that k_{z, R} ≠ 0 regardless of the value of ω. This apparent contradiction between the standing and propagating cases can be understood as follows. An oscillatory standing wave can be interpreted as the superposition of two propagating waves with the same frequency running in opposite directions. However, such a representation does not work for a perturbation fixed in space and evanescent in time. The presence of static evanescent perturbations is a peculiar result only obtained when k_z is real and ω is purely imaginary, a situation that cannot be described with propagating waves. For this reason evanescent solutions are absent from the present study of propagating waves.

3.2.1. Approximate Expressions and Comparison with Numerical Results

Going back to Equations (35) and (36), we realize that it is possible to obtain simpler expressions of $k_{z,\rm R}^2$ and $k_{z,\rm I}^2$ by taking advantage of the fact that the second term within the square root of Equation (35) is always much smaller than unity. We can approximate Equations (35) and (36) as

$$\begin{eqnarray} k_{z,\rm R}^2 &\approx & \frac{\omega ^2}{c_{\mathrm{A}}^2} \frac{\omega ^2 + (1+\chi)\nu _{\mathrm{ni}}^2}{\omega ^2 + \nu _{\mathrm{ni}}^2} \nonumber \\ &&+\, \frac{1}{4} \frac{\omega ^2}{c_{\mathrm{A}}^2} \frac{\chi ^2\nu _{\mathrm{ni}}^2\omega ^2}{\left(\omega ^2 + \nu _{\mathrm{ni}}^2\right)\left(\omega ^2 + (1+\chi)\nu _{\mathrm{ni}}^2\right)}, \end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray} k_{z,\rm I}^2 &\approx & \frac{1}{4} \frac{\omega ^2}{c_{\mathrm{A}}^2} \frac{\chi ^2\nu _{\mathrm{ni}}^2\omega ^2}{\left(\omega ^2 + \nu _{\mathrm{ni}}^2\right)\left(\omega ^2 + (1+\chi)\nu _{\mathrm{ni}}^2\right)}. \end{eqnarray} \tag{ 38 }$$

As for standing waves we evaluate the analytic solutions (Equations (37) and (38)) in the various limits of ν_ni. When ν_ni ≪ ω, Equations (37) and (38) simplify to

$$\begin{eqnarray} k_{z,\rm R} &\approx & {\pm }\sqrt{ \frac{\omega ^2}{c_{\mathrm{A}}^2} + \frac{\chi ^2\nu _{\mathrm{ni}}^2}{4c_{\mathrm{A}}^2}} \approx \pm \frac{\omega }{c_{\mathrm{A}}}, \end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray} k_{z,\rm I} &\approx & {\pm }\frac{\chi \nu _{\mathrm{ni}}}{2c_{\mathrm{A}}}.\qquad\qquad\quad\;\; \end{eqnarray} \tag{ 40 }$$

Hence, we find that k_{z, R} takes the same value as in the fully ionized case, and k_{z, I} is independent of ω. Conversely when ν_ni ≫ ω we find

$$\begin{eqnarray} k_{z,\rm R} &\approx & {\pm }\sqrt{ \frac{\omega ^2}{c_{\mathrm{A}}^2}(1+\chi) + \frac{1}{4}\frac{\omega ^2}{c_{\mathrm{A}}^2}\frac{\chi ^2\omega ^2}{(1+\chi)^2\nu _{\mathrm{in}}^2} } \nonumber \\ &\approx & {\pm }\frac{\omega }{c_{\mathrm{A}}}\sqrt{1+\chi }, \end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray} k_{z,\rm I} &\approx & {\pm }\frac{\chi \omega ^2}{2(1+\chi)c_{\mathrm{A}}\nu _{\mathrm{in}}}. \qquad\qquad\quad\, \end{eqnarray} \tag{ 42 }$$

We find the presence of the factor $\sqrt{1+\chi }$ in the approximate expression of k_{z, R}, while now k_{z, I} is proportional to ω². Hence, high-frequency waves are more efficiently damped than low-frequency waves.

Now we solve the dispersion relation (Equation (14)) numerically and compare the numerical solutions with the analytic approximations of k_{z, R} and k_{z, I} (Equations (37) and (38)). These results are shown in Figure 4 for a particular set of parameters given in the caption of the figure. As in the standing case, the agreement between numerical and analytic results is very good. We have reapplied these computations using various values of χ (not shown here for simplicity) and the analytic results are always in accordance with the numerical ones.

Zoom In Zoom Out Reset image size

3.2.2. Phase Velocity and Group Velocity

Here we compute the phase and group velocities of the propagating Alfvén waves. The phase velocity, v_ph, is defined as

$$\begin{equation} {\bf v}_{\rm ph} = \frac{\omega }{k} \hat{e}_k, \end{equation} \tag{ 43 }$$

with $\hat{e}_k$ denoting the unit vector in the direction of the wave vector. From the dispersion relation (Equation (14)) we get

$$\begin{equation} \frac{\omega }{k} = c_{\mathrm{A}}\cos \theta \sqrt{\frac{\omega +i\nu _{\mathrm{ni}}}{\omega +i(1+\chi)\nu _{\mathrm{ni}}}}. \end{equation} \tag{ 44 }$$

The phase velocity is a complex quantity. Its real part is related to the propagation speed of the wave, while its imaginary part is related to the wave attenuation. Both real and imaginary parts of v_ph have the same dependence on θ, meaning that the damping rate is independent of the direction of propagation. Ion–neutral collisions do not affect the dependence of v_ph on θ. From Equation (44) we can define the effective Alfvén velocity, $\tilde{c}_{\rm A}$, as

$$\begin{equation} \tilde{c}_{\rm A} = c_{\mathrm{A}}\sqrt{\frac{\omega +i\nu _{\mathrm{ni}}}{\omega +i(1+\chi)\nu _{\mathrm{ni}}}}. \end{equation} \tag{ 45 }$$

The group velocity, v_gr, is the propagation velocity of a wave packet and, therefore, of the wave energy. It is defined as

$$\begin{equation} {\bf v}_{\rm gr} = \nabla _{\bf k} \omega = \frac{\partial \omega }{\partial k_x} \hat{e}_x + \frac{\partial \omega }{\partial k_y} \hat{e}_y + \frac{\partial \omega }{\partial k_z} \hat{e}_z. \end{equation} \tag{ 46 }$$

Using the dispersion relation (Equation (14)) the group velocity can be more easily computed as

$$\begin{equation} {\bf v}_{\rm gr} = - \frac{1}{\partial \mathcal {D} / \partial \omega } \left(\frac{\partial \mathcal {D}}{\partial k_x} \hat{e}_x + \frac{\partial \mathcal {D}}{\partial k_y} \hat{e}_y + \frac{\partial \mathcal {D}}{\partial k_z} \hat{e}_z \right) \end{equation} \tag{ 47 }$$

that gives

$$\begin{eqnarray} {\bf v}_{\rm gr} &=& \frac{2 k_z c_{\mathrm{A}}^2 (\omega + i\nu _{\mathrm{ni}})}{3\omega ^2 - k_z^2c_{\mathrm{A}}^2 + i2(1+\chi)\nu _{\mathrm{ni}}\omega } \hat{e}_z \nonumber \\ &=& 2 c_{\mathrm{A}}\frac{\left(\omega +i\nu _{\mathrm{ni}}\right)^{3/2}\left(\omega + i(1+\chi)\nu _{\mathrm{ni}}\right)^{1/2}}{2\omega ^2-2(1+\chi)\nu _{\mathrm{ni}}^2+i(1+\chi)\nu _{\mathrm{ni}}\omega } \hat{e}_z. \end{eqnarray} \tag{ 48 }$$

As it happens for the phase velocity, the group velocity is also complex. Since the damping of propagating Alfvén waves depends on ω, the components of the wave packet with higher ω are more damped than the components with low ω. Hence, the concept of the group velocity as the velocity at which the whole wave packet propagates becomes obsolete when there is dissipation. However, following Muschietti & Dum (1993) it is still possible to define an effective group velocity as the velocity at which the center of the wave packet propagates. This effective group velocity is a time-dependent combination of the real and imaginary parts of v_gr and depends on the particular form of the wave packet (see extensive details in Muschietti & Dum 1993).

In the limits of low and high collision frequency the damping is weak, so that the imaginary part of v_gr can be neglected compared to the real part. When ν_ni ≪ ω, we find that ${\bf v}_{\rm gr} \approx c_{\mathrm{A}}\hat{e}_z$ as in the ideal case, while when ν_ni ≫ ω, the group velocity is ${\bf v}_{\rm gr} \approx c_{\mathrm{A}}/ \sqrt{1+\chi }\hat{e}_z$. Note again the presence of the factor $\sqrt{1+\chi }$.

In Section 3, we have followed a normal mode approach and have assumed that the temporal dependence of the perturbations is of the form exp (− iωt). Here we go beyond the normal mode analysis and look for general time-dependent solutions to Equations (9) and (10). To this end we follow two different approaches. First, we solve the time-dependent problem analytically by means of the Laplace transform. Later we use the numerical code MolMHD (Bona et al. 2009) to numerically evolve the full set of basic Equations (1)–(5). Finally, these two independent results are compared. In this section we restrict ourselves to standing waves and assume a real and positive k_z. We refer the reader to Roberge & Ciolek (2007) for the problem of driven propagating waves in a two-fluid plasma.

4.1. Analytic Solution

Here we look for general time-dependent solutions to the coupled Equations (9) and (10) using the Laplace transform. We compute the Laplace transform in time of the vorticity perturbations Γ_i and Γ_n as

$$\begin{equation} \tilde{\Gamma }_{\beta }(s) = \mathcal {L}[ \Gamma _\beta ] = \int _0^\infty \Gamma _\beta (t) \exp \left(-s t \right) {\rm d}t, \end{equation} \tag{ 49 }$$

where s is the transformed variable and β denotes either "i" or "n." The Laplace transforms of the temporal derivatives of vorticity are given by

$$\begin{eqnarray} \mathcal {L}\left[ \frac{\partial \Gamma _\beta }{\partial t} \right] &=& s \tilde{\Gamma }_{\beta }(s) - \Gamma _{0,\beta }, \end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray} \mathcal {L}\left[ \frac{\partial ^2 \Gamma _\beta }{\partial t^2} \right] &=& s^2 \tilde{\Gamma }_{\beta }(s) - s \Gamma _{0,\beta } - \Gamma ^{\prime }_{0,\beta }, \end{eqnarray} \tag{ 51 }$$

where Γ_{0, β} and $\Gamma ^{\prime }_{0,\beta }$ denote the value of Γ_β and its temporal derivative at t = 0, respectively. From Equations (9) and (10), we deduce that the temporal derivatives of the vorticity perturbations at t = 0 satisfy

$$\begin{eqnarray} \Gamma ^{\prime }_{0,\rm i} &=& {-}\chi \nu _{\mathrm{ni}}(\Gamma _{0,\rm i} - \Gamma _{0,\rm n} ), \end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray} \Gamma ^{\prime }_{0,\rm n} &=& \nu _{\mathrm{ni}}(\Gamma _{0,\rm i} - \Gamma _{0,\rm n}). \end{eqnarray} \tag{ 53 }$$

At the present stage, we leave Γ_{0, i} and Γ_{0, n} unspecified.

We apply the Laplace transform to Equations (9) and (10). From the transformed Equation (10) we can express $\tilde{\Gamma }_{\rm n}$ in terms of $\tilde{\Gamma }_{\rm i}$ as

$$\begin{equation} \tilde{\Gamma }_{\rm n}(s) = \frac{1}{s + \nu _{\mathrm{ni}}} \left(\nu _{\mathrm{ni}}\tilde{\Gamma }_{\rm i}(s) + \Gamma _{0, \rm n} \right). \end{equation} \tag{ 54 }$$

We insert this expression in the transformed Equation (9) and obtain the expression for $\tilde{\Gamma }_{\rm i}$ as

$$\begin{equation} \tilde{\Gamma }_{\rm i}(s) = \frac{\mathcal {B}(s)}{\mathcal {D}^*(s)}, \end{equation} \tag{ 55 }$$

where the functions $\mathcal {B}(s)$ and $\mathcal {D}^*(s)$ are defined as

$$\begin{eqnarray} \mathcal {B}(s) &=& s [(s + \nu _{\mathrm{ni}})\Gamma _{0, \rm i} + \chi \nu _{\mathrm{ni}}\Gamma _{0, \rm n} ], \end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray} \mathcal {D}^*(s) &=& s^3 + (1+\chi ) \nu _{\mathrm{ni}}s^2 + k_z^2c_{\mathrm{A}}^2 s + \nu _{\mathrm{ni}}k_z^2c_{\mathrm{A}}^2. \end{eqnarray} \tag{ 57 }$$

The function $\mathcal {B}(s)$ depends on the initial conditions and so it contains information about how the waves are excited. Conversely, the function $\mathcal {D}^*(s)$ is independent of the initial conditions. The function $\mathcal {D}^*(s)$ plays a very important role because it corresponds to the dispersion function. It tells us how the plasma behaves after the excitation has taken place. We note that the expression of $\mathcal {D}^*(s)$ (Equation (57)) coincides with the left-hand side of Equation (16), which was obtained from the normal mode dispersion relation (Equation (14)) after performing the change of variable ω = −is. Hence the normal modes are consistently recovered from the present Laplace transform approach by setting $\mathcal {D}^*(s) = 0$.

To compute the vorticity perturbation, Γ_i, in the actual temporal domain we perform the inverse Laplace transform of Equation (55), namely

$$\begin{equation} \Gamma _{\rm i}(t) = \mathcal {L}^{-1} \left[ \tilde{\Gamma }_{\rm i}(s) \right] = \mathcal {L}^{-1} \left[ \frac{\mathcal {B}(s)}{\mathcal {D}^*(s)} \right]. \end{equation} \tag{ 58 }$$

The inverse Laplace transform of Equation (58) can be evaluated by using the technique of partial fraction decomposition (see, e.g., Dyke 1999) if the roots of the equation $\mathcal {D}^*(s) = 0$ are known. The zeros of $\mathcal {D}^*(s)$ are precisely the normal modes studied in Section 3, so that we can take advantage of the approximate expressions found in Section 3. Once Γ_i(t) is known, the corresponding expression of Γ_n(t) can be obtained using Equation (54) and applying the convolution theorem, namely

$$\begin{equation} \Gamma _{\rm n}(t) = \Gamma _{0, \rm n} e^{-\nu _{\mathrm{ni}}t} + \nu _{\mathrm{ni}}\int _0^t \Gamma _{\rm i}(\tau) \exp \left[-\nu _{\mathrm{ni}}\left(t-\tau \right) \right] {\rm d}\tau. \end{equation} \tag{ 59 }$$

From Equations (58) and (59) we obtain Γ_i(t) and Γ_n(t), namely

$$\begin{eqnarray} \Gamma _{\rm i}(t) &=& A_1 \exp \left(\epsilon t \right) \nonumber \\ &&+ \left[ A_2 \cos \left(\omega _{\rm R} t\right) + A_3 \sin \left(\omega _{\rm R} t\right) \right]\exp \left(\omega _{\rm I} t \right), \end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray} \Gamma _{\rm n}(t) &=& C_1 \exp \left(\epsilon t \right) + C_2 \exp \left(-\nu _{\mathrm{ni}}t \right) \nonumber \\ && + \left[ C_3 \cos \left(\omega _{\rm R} t\right) + C_4 \sin \left(\omega _{\rm R} t\right) \right]\exp \left(\omega _{\rm I} t \right), \end{eqnarray} \tag{ 61 }$$

with the approximate expressions of ω_R, ω_I, and given in Equations (21)–(23) and the coefficients A₁–A₃ and C₁–C₄ given in the Appendix. In the expression of ω_R (Equation (21)) the + sign has to be taken. We must recall that the expressions of ω_R, ω_I, and given in Equations (21)–(23) are computed in the weak damping approximation.

In this analysis we have worked with the vorticity perturbations. However, we can also write the solution to the initial-value problem using velocity perturbations. To do so we move, for convenience, to a reference frame in which k_y = 0. Thus, the motions of Alfvén waves are polarized in the y-direction and from Equations (7) and (8) we get that the vorticity perturbations are proportional to the y-component of velocities, namely v_{i, y} and v_{n, y}. Therefore, the expressions for the temporal evolution of v_{i, y} and v_{n, y} are the same as those given in Equations (60) and (61) for Γ_i and Γ_n, respectively. We just need to replace Γ_{0, i} by v_{i, 0} and Γ_{0, n} by v_{n, 0}, where v_{i, 0} and v_{n, 0} denote the values of the ion and neutral velocities at t = 0, respectively.

4.1.1. Solution in the Absence of Magnetic Field

We start the study of some interesting cases by considering the situation in which there is no magnetic field. We set B = 0, and therefore c_A = 0. This situation was studied by Vranjes et al. (2008). Equations (60) and (61) become

$$\begin{eqnarray} \Gamma _{\rm i}(t) &=& \frac{\Gamma _{0, \rm i} + \chi \Gamma _{0, \rm n}}{1+\chi } \nonumber \\ &&+ \frac{\chi \left(\Gamma _{0, \rm i} - \Gamma _{0, \rm n} \right) }{1+\chi }\exp \left[ - \left(1 + \chi \right) \nu _{\mathrm{ni}}t \right], \end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray} \Gamma _{\rm n}(t) &=& \frac{\Gamma _{0, \rm i} + \chi \Gamma _{0, \rm n}}{1+\chi } \nonumber \\ &&- \frac{\left(\Gamma _{0, \rm i} - \Gamma _{0, \rm n} \right) }{1+\chi }\exp \left[ - \left(1 + \chi \right) \nu _{\mathrm{ni}}t \right]. \end{eqnarray} \tag{ 63 }$$

These solutions agree with those found by Vranjes et al. (2008). After an initial phase, i.e., when (1 + χ)ν_nit ≫ 1, both Γ_i and Γ_n relax to the same constant value, which corresponds to a weighted average of the initial conditions, namely

$$\begin{equation} \Gamma _{\rm i}(t) \approx \Gamma _{\rm n}(t) \approx \frac{\Gamma _{0, \rm i} + \chi \Gamma _{0, \rm n}}{1+\chi } = \frac{\rho _{\rm i}\Gamma _{0, \rm i} + \rho _{\rm n}\Gamma _{0, \rm n}}{\rho _{\rm i}+\rho _{\rm n}}. \end{equation} \tag{ 64 }$$

In the absence of magnetic field no oscillatory solutions are found. Vorticity perturbations remain constant after the relaxation phase.

4.1.2. Low Collision Frequency

Here we incorporate again the magnetic field but we consider the case of low collision frequency compared to the oscillation frequency, i.e., ν_ni ≪ ω_R. Equations (60) and (61) simplify to

$$\begin{eqnarray} \Gamma _{\rm i}(t) & \approx & \Gamma _{0, \rm i} \cos \left(k_z c_{\mathrm{A}}t \right) \exp \left(- \frac{\chi \nu _{\mathrm{ni}}}{2} t \right), \end{eqnarray} \tag{ 65 }$$

$$\begin{eqnarray} \Gamma _{\rm n}(t) & \approx & \Gamma _{0, \rm n} \exp \left(- \nu _{\mathrm{ni}}t \right). \end{eqnarray} \tag{ 66 }$$

When ν_ni ≪ ω_R the vorticity perturbations of the ionized fluid and the neutral fluid are essentially decoupled from each other. After the excitation, vortical motions in the ionized fluid oscillate at the Alfvén frequency, k_zc_A, and with a damping rate proportional to ν_ni. These are weakly damped Alfvén waves. Conversely, vortical disturbances in the neutral fluid are evanescent in time. There is no driving force for vortical motions in the neutral fluid, and so there is no oscillatory behavior.

4.1.3. High Collision Frequency

Now we turn to the limit of high collision frequency compared to the oscillation frequency. This is the realistic case for many astrophysical applications and deserves special attention. We compute the limit of Equations (60) and (61) when ν_ni ≫ ω_R. We obtain

$$\begin{eqnarray} \Gamma _{\rm i}(t) & \approx & \frac{ \Gamma _{0, \rm i} + \chi \Gamma _{0, \rm n}}{1+\chi } \cos \left(\frac{k_zc_{\mathrm{A}}}{\sqrt{1+\chi }} t \right) \nonumber \\ &&\times \exp \left(-\frac{\chi }{2\left(1+\chi \right)^2 } \frac{k_z^2 c_{\rm A}^2 }{\nu _{\mathrm{ni}}} t \right) \nonumber \\ &&+ \frac{\chi \left(\Gamma _{0, \rm i} - \Gamma _{0, \rm n} \right) }{1+\chi }\exp \left[ - \left(1 + \chi \right) \nu _{\mathrm{in}}t \right], \end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray} \Gamma _{\rm n}(t) & \approx & \frac{ \Gamma _{0, \rm i} + \chi \Gamma _{0, \rm n}}{1+\chi }\cos \left(\frac{k_zc_{\mathrm{A}}}{\sqrt{1+\chi }} t \right) \nonumber \\ &&\times \exp \left(-\frac{\chi }{2\left(1+\chi \right)^2 } \frac{k_z^2 c_{\rm A}^2 }{\nu _{\mathrm{ni}}} t \right) \nonumber \\ &&- \frac{\Gamma _{0, \rm i} - \Gamma _{0, \rm n} }{1+\chi }\exp \left[ - \left(1 + \chi \right) \nu _{\mathrm{in}}t \right]. \end{eqnarray} \tag{ 68 }$$

By comparing Equations (67) and (68) with the equivalent solutions in the absence of magnetic field (Equations (62) and (63)) we find that in both cases there is a relaxation phase whose timescale is determined by the collision frequency. In the absence of magnetic field (Equations (62) and (63)), the vorticity perturbations remain constant after the relaxation phase, while in the presence of magnetic field (Equations (67)–(68)) the perturbations of both ions and neutrals oscillate together as a single fluid at the modified Alfvén frequency $k_z c_{\mathrm{A}}/\sqrt{1+\chi }$. The oscillations are exponentially damped in time.

We define Δ(t) as the difference of Γ_i(t) and Γ_n(t) computed from Equations (67) and (68), namely

$$\begin{equation} \Delta (t) \equiv \Gamma _{\rm i}(t) - \Gamma _{\rm n}(t) \approx \left(\Gamma _{0, \rm i} -\Gamma _{0, \rm n} \right) \exp \left(- \frac{t}{\tau _{\rm rel}}\right), \end{equation} \tag{ 69 }$$

with τ_rel the relaxation timescale given by

$$\begin{equation} \tau _{\rm rel} = \frac{1}{\left(1+\chi \right)\nu _{\mathrm{ni}}} = \frac{1}{\nu _{\mathrm{in}}+ \nu _{\mathrm{ni}}}. \end{equation} \tag{ 70 }$$

The function Δ(t) informs us about the relaxation phase. The function Δ(t) → 0 when t ≫ τ_rel. Importantly, Δ(t) = 0 when $\Gamma _{0, \rm i} =\Gamma _{0, \rm n}$, i.e., there is no relaxation phase when the initial disturbance perturbs ions and neutrals in the same way. In such a case, both ions and neutrals oscillate together with amplitude equal to the initial condition. When $\Gamma _{0, \rm i} \ne \Gamma _{0, \rm n}$ the amplitude of the Alfvénic oscillations is the weighted average of the initial conditions (Equation (64)). The weights correspond to the densities of each fluid.

4.2. Numerical Solution

Here we use the numerical code MolMHD to evolve in time the basic Equations (1)–(5) in a fully numerical way. The numerical code (see Bona et al. 2009, for details about the scheme) uses the method of lines for the discretization of the variables, and the time and space variables are treated separately. For the temporal part, a fourth-order Runge–Kutta method is used, whereas for the space discretization a finite difference scheme with a fourth-order centered stencil is used. For a given spatial resolution, the time step is selected so as to satisfy the Courant condition.

The MolMHD code evolves in time the plasma and magnetic field perturbations after an initial condition. In the present application the perturbations are assumed invariant in the x- and y-directions, so that perturbations only depend on the z-direction and Alfvén waves are strictly polarized in the y-direction. The only non-zero perturbations are the components of velocity and magnetic field in the y-direction. The numerical integration of Equations (1)–(5) is done in the interval z ∈ [ − 10L, 10L], where L is an arbitrary length scale. We use a uniform grid with 1001 grid points. Since we study standing waves, the boundary conditions used at z = ±10L are that the velocity perturbations are fixed to zero. The boundary conditions for the other wave perturbations are that their spatial derivatives are set to zero. The initial condition at t = 0 for the velocity of ions, v_{i, y}, is

$$\begin{equation} v_{{\rm i},y}= v_{{\rm i},0} \cos \left(\frac{\pi }{20L}z \right), \end{equation} \tag{ 71 }$$

which corresponds to the fundamental standing mode in our domain with dimensionless wavenumber k_zL = π/20. In the fully ionized case, the dimensionless frequency of the fundamental mode is ωL/c_A = π/20 ≈ 0.157. Since we are dealing with linear perturbations, we express the velocity perturbation in arbitrary units and set v_{i, 0} = 1.

Figure 5 shows the velocity perturbation of ions and neutrals at z = 0 computed numerically with the MolMHD code for various values of the neutral–ion collision frequency. We compare the numerical results with the analytic ones (Equations (60) and (61)). The results in the top row of Figure 5 (panels (a)–(c)) are obtained when the velocity of neutrals at t = 0 is the same as that of ions (Equation (71)), while in the bottom row of Figure 5 (panels (d)–(f)) neutrals are initially at rest. In the following paragraphs we discuss the results displayed in Figure 5 and relate the time-dependent solutions with the normal modes investigated in Section 3.

• 1.  
  In Figures 5(a) and (d), we use ν_niL/c_A = 0.01 so that the collision frequency is an order of magnitude lower than the ideal Alfvén frequency. Ions and neutrals behave almost independently. Regardless of the initial condition for neutrals, ions oscillate at the Alfvén frequency k_zc_A and are damped due to collisions. The oscillations of the ionized fluid act as a periodic forcing on neutrals, although the oscillations generated in the neutral fluid are of much lower amplitude compared to those in the ionized fluid. The dominant behavior in neutrals when the initial perturbation is non-zero is the exponential decay of the initial perturbation. This behavior is consistent with Equations (65) and (66). In relation with the normal modes of Section 3, we see that when the collision frequency is low the oscillatory mode is mostly excited in the ionized fluid, while the evanescent mode dominates the neutral fluid behavior. Hence, we can understand the physical reason for the existence of the evanescent mode as the way in which vorticity perturbations in the neutral fluid decay in time due to collisions (see also Zaqarashvili et al. 2011b).
• 2.  
  In Figures 5(b) and (e), we increase the collision frequency to ν_niL/c_A = 0.1. The collision frequency and the ideal Alfvén frequency are of the same order of magnitude. Now both ions and neutrals display strongly damped oscillations and it is not possible to relate the behavior of each fluid with one particular normal mode. Both oscillatory and evanescent normal modes are excited in both fluids and the observed behavior is the result of the joint effect of both oscillatory and evanescent modes. We note that in Figures 5(b) and (e) the analytic solution of the time-dependent problem does not exactly follow the full numerical result. The source of the discrepancy is in the approximate expressions of ω_R, ω_I, and (Equations (21)–(23)) used in the analytic solution of the initial-value problem (Equations (60) and (61)). These approximate expressions were obtained in the case of weak damping, which obviously does not apply to the situation of Figures 5(b) and (e). However, it is remarkable that the analytic solution still captures the overall behavior correctly even beyond the weak damping limit. It is worth noting that instead of using the approximate values of ω_R, ω_I, and , we could use their exact values obtained by numerically solving the dispersion relation. In such a case the agreement between analytic and numerical results is excellent (not shown here for simplicity).
• 3.  
  Finally, we use ν_niL/c_A = 1 in Figures 5(c) and (f) so that the collision frequency is an order of magnitude higher than the ideal Alfvén frequency. We find that ions and neutrals behave as a single fluid and oscillate together at the modified Alfvén frequency $k_z c_{\mathrm{A}}/\sqrt{1+\chi }$. Again, both numerical and analytic results are in excellent agreement. The common amplitude of the y-components of velocity of ions and neutrals, namely $\hat{v}_{y}$, is

  $$\begin{equation} \hat{v}_{y} = \frac{\rho _{\rm i}v_{{\rm i},0} + \rho _{\rm n}v_{{\rm n},0}}{\rho _{\rm i}+ \rho _{\rm n}}. \end{equation} \tag{ 72 }$$

  The oscillations are weakly damped. When the initial velocity of neutrals is the same as that of ions (Figure 5(c)), both fluids oscillate together from t = 0. The evanescent mode is not excited when v_{i, 0} = v_{n, 0}. Conversely, when neutrals are initially at rest (Figure 5(f)), the single-fluid behavior begins after a short relaxation phase. Although very brief, the relaxation phase can be seen in Figure 5(f) near t = 0 as a sudden decrease of v_{i, y} and increase of v_{n, y}. Figure 6 shows the same results displayed in Figure 5(f) but near t = 0 in order to observe the relaxation phase in detail. The behavior of the velocity perturbations is governed by the evanescent normal mode during the relaxation phase and by the oscillatory normal mode afterward. This is consistent with Equations (67) and (68). In this case and since v_{n, 0} = 0, the amplitude of the joint oscillations after the relaxation phase (Equation (72)) is

  $$\begin{equation} \hat{v}_{y} = \frac{\rho _{\rm i}}{\rho _{\rm i}+ \rho _{\rm n}} v_{{\rm i},0} = \frac{1 }{1 + \chi } v_{{\rm i},0}. \end{equation} \tag{ 73 }$$

  Hence, the larger the ionization ratio χ, the lower the amplitude.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this section we perform an application of the previous theoretical analysis to a particular astrophysical plasma, namely the solar atmosphere. We specialize in solar atmospheric plasma because of two main reasons. The plasma in the coolest parts of the solar atmosphere, i.e., the photosphere and the low chromosphere, is weakly ionized. In addition, recent observations have shown the ubiquitous presence of Alfvénic waves in the solar atmosphere (e.g., De Pontieu et al. 2007; Tomczyk et al. 2007; McIntosh et al. 2011; Okamoto & De Pontieu 2011; De Pontieu et al. 2012). It is therefore interesting to apply the theory developed in the present paper to the solar case. Many works have focused on the damping of chromospheric waves due to ion–neutral collisions and its importance for plasma heating (see, e.g., Haerendel 1992; De Pontieu et al. 2001; Khodachenko et al. 2004; Leake et al. 2005; Soler et al. 2012; Zaqarashvili et al. 2013, among others). In the present application we do not discuss damping but investigate other effects on Alfvén waves caused by partial ionization, namely the presence of cutoff wavelengths of standing waves and the modification of the effective Alfvén velocity. For results about wave damping, the reader is referred to the papers cited above.

We consider a simplified model for the solar atmospheric plasma. The variation of physical parameters with height from the photosphere to the base of the corona is taken from the quiet-Sun model C of Vernazza et al. (1981), hereafter VALC model. To compute the friction coefficient, α_in, we ignore the influence of heavier species and consider only hydrogen. Hence, the expression of α_in is given in Equation (6), where we take σ_in ≈ 10⁻²⁰ m² according to the estimation by Zaqarashvili et al. (2013) for the case of direct elastic collisions (see also Braginskii 1965). Figure 7 shows the inverse of the relaxation time of the ion–neutral coupling, $\tau _{\rm rel}^{-1}$ (Equation (72)), as a function of height in the low solar atmosphere. The very large values of $\tau _{\rm rel}^{-1}$, i.e., very short τ_rel, point out that in the solar atmosphere ions and neutrals are very efficiently coupled and, for practical purposes, they can be considered as a single fluid. The single-fluid approach is usually followed in studies focused on wave damping (see the references in the previous paragraph). For comparison, we also plot in Figure 7 the ion–neutral, ν_in, and neutral–ion, ν_ni, collision frequencies (see also similar plots in De Pontieu et al. 2001). We find that the value of τ_rel is mainly determined by ν_in. Note that in Figure 7 the solid and dotted lines are superimposed except at large heights.

Zoom In Zoom Out Reset image size

Concerning the magnetic field strength, we consider the model by Leake & Arber (2006) of a vertical chromospheric magnetic flux tube expanding with height, whose form is given by

$$\begin{equation} B = B_{\rm ph} \left(\frac{\rho }{\rho _{\rm ph}} \right)^\beta, \end{equation} \tag{ 74 }$$

where ρ = ρ_i + ρ_n is the total density, B_ph and ρ_ph are the magnetic field strength and the total density, respectively, at the photospheric level, and β = 0.3 is an empirical exponent. We use ρ_ph = 2.74 × 10⁻⁴ kg m⁻³ from the VALC model and B_ph = 1.5 kG. The dependence of ρ with height is determined by the VALC model. With this choice of parameters the magnetic field strength decreases with height so that B ≈ 100 G at 1000 km and B ≈ 20 G at 2000 km above the photosphere.

The physical parameters in the simplified model of the low solar atmosphere used here depend on the vertical direction, z, while in the theoretical analysis of the previous sections all the parameters are taken constant in z. However, it is possible to apply the expressions derived before for a homogeneous plasma to the present stratified case if the wavelength in the z-direction, λ_z = 2π/k_z, is much shorter than the length scale of the variation of the effective Alfvén velocity. If such a restriction is fulfilled, we can perform a local analysis and use the physical parameters at a given height in the expressions derived for homogeneous plasma. We follow this approach.

Another assumption made in this section is that the amplitudes of the waves are small enough for the linear analysis to remain valid. This is a reasonable assumption for chromospheric waves. A parameter that can quantify nonlinearity of the waves is the ratio of the wave velocity amplitude to the local Alfvén velocity. The median of the velocity amplitudes of the chromospheric waves detected by Okamoto & De Pontieu (2011) is 7.4 km s⁻¹. This value is, at least, an order of magnitude lower than the expected value of the Alfvén velocity in the chromosphere.

We start by computing the cutoff region of wavenumbers, k_z, of standing waves (Equation (20)) as a function of height. Instead of k_z, we show in Figure 8(a) the corresponding wavelength, λ_z. The shaded area in Figure 8(a) indicates the range of wavelengths where oscillatory standing modes are not possible. At low heights above the photosphere the cutoff wavelengths are very short. The cutoff wavelengths increase with height and become of the order of a few kilometers at heights between 500 km and 1500 km, approximately. Then the cutoff region disappears at a height of 1600 km above the photosphere, approximately, because the ionization ratio, χ, decreases with height as the plasma becomes more and more ionized and the threshold value χ = 8 is reached at that height.

Zoom In Zoom Out Reset image size

Next we turn to the computation of the effective Alfvén velocity (Equation (45)). Since this quantity depends on the wave frequency we use a frequency of 22 mHz, which corresponds to the dominant frequency of the chromospheric Alfvénic waves detected by Okamoto & De Pontieu (2011). We display in Figure 8(b) the real part of the effective Alfvén velocity as a function of height. In Figure 8(b) we also show the corresponding Alfvén velocity computed using its classic definition that depends on the ion density only. We find that in the low chromosphere neutral–ion collisions are crucial to decrease the effective Alfvén velocity between one and two orders of magnitude with respect to the expected value taking into account the ion density only.

Now we can check the assumption that λ_z is much shorter than the length scale of the variation of the effective Alfvén velocity. On the one hand from Figure 8(b), we get that the effective Alfvén velocity varies about an order of magnitude in 2000 km, approximately. On the other hand, the largest cutoff wavelengths displayed in Figure 8(a) are of the order of a few tens of kilometers. Hence, the use of a local analysis in the present application is justified.

In the solar context, the dissipation of Alfvénic waves may play an important role for the heating of the atmospheric plasma (see, e.g., De Pontieu et al. 2007; McIntosh et al. 2011). The implications of partial ionization for the calculations of energy carried by Alfvén waves were discussed by, e.g., Vranjes et al. (2008) and Tsap et al. (2011). Both papers give expressions for the energy flux of Alfvén waves in partially ionized plasmas. However, the equations provided by Vranjes et al. (2008) and Tsap et al. (2011) are different. The expression of Vranjes et al. (2008) includes the factor (ρ_i/ρ_n)², which is absent from the expression of Tsap et al. (2011). Hence, the conclusions of these two papers regarding the impact of partial ionization are in apparent contradiction. On the one hand, Vranjes et al. (2008) explained that due to the factor (ρ_i/ρ_n)² the energy flux in weakly ionized plasmas becomes orders of magnitude smaller compared to the fully ionized case. On the other hand, since Tsap et al. (2011) lacked that factor they argued that the energy flux is independent of the ionization degree and obtained the same expression as for fully ionized plasmas. Here we shall try to solve this apparent contradiction between the results of Vranjes et al. (2008) and Tsap et al. (2011).

For application to the solar atmosphere we take the limit of high collision frequency. We choose a reference frame in which k_y = 0. After the relaxation phase, i.e., for t ≫ τ_rel, and neglecting the weak damping due to ion–neutral collisions, the y-components of velocity of ions and neutrals oscillate together with frequency, $k_z c_{\mathrm{A}}/\sqrt{1+\chi }$, and amplitude, $\hat{v}_{y}$, given in Equation (72). The energy flux along magnetic field lines, S_z, calculated using the Poynting vector (see, e.g., Walker 2005) is

$$\begin{equation} S_z= \frac{1}{2 \mu } \hat{E}_{x} \hat{b}_{y}, \end{equation} \tag{ 75 }$$

where $\hat{E}_{x}$ and $\hat{b}_{y}$ denote the amplitudes of the x-component of the electric field and the y-component of the magnetic field perturbation. The electric field is

$$\begin{equation} {\bf E} = - {\bf v}_{\rm i} \times {\bf B}, \end{equation} \tag{ 76 }$$

so that $\hat{E}_{x}\sim \hat{v}_{y} B$. On the other hand, the magnetic field perturbation is governed by Equation (3), from where we get $\hat{b}_{y}\sim \hat{v}_{y} B \sqrt{1+\chi }/ c_{\mathrm{A}}$. Hence, the energy flux is

$$\begin{equation} S_z = \frac{\sqrt{1+\chi }}{c_{\mathrm{A}}}\frac{B^2}{2\mu } \hat{v}_{y}^2= \frac{B}{2\sqrt{\mu }} \sqrt{\rho } \, \hat{v}_{y}^2, \end{equation} \tag{ 77 }$$

with ρ = ρ_i + ρ_n the total plasma density. Consistently we find the classical expression of the energy flux of Alfvén waves (see, e.g., Walker 2005), with the only difference that the total plasma density replaces ion density and that the velocity amplitude after the relaxation phase, $\hat{v}_{y}$, is used. Since the relaxation phase is extremely short for realistic collision frequencies (Figure 7), $\hat{v}_{y}$ is the velocity amplitude that an observer would actually measure. In a fully ionized plasma, $\hat{v}_{y} = v_{{\rm i},0}$ so that

$$\begin{equation} S_z = \frac{B}{2\sqrt{\mu }} \sqrt{\rho } \, v_{{\rm i},0}^2. \end{equation} \tag{ 78 }$$

Up to here there is no discrepancy between the analysis of Vranjes et al. (2008) and Tsap et al. (2011). The discrepancy arises when the energy flux is written in terms of the initial velocity amplitudes of ions and neutrals. We insert Equation (72) in Equation (77) and rewrite the energy flux in terms of the initial velocity amplitudes, namely

$$\begin{equation} S_z = \frac{B}{2\sqrt{\mu }} \sqrt{\rho } \left(\frac{\rho _{\rm i}v_{{\rm i},0} + \rho _{\rm n}v_{{\rm n},0}}{\rho _{\rm i}+ \rho _{\rm n}} \right)^2. \end{equation} \tag{ 79 }$$

In weakly ionized plasmas ρ_n ≫ ρ_i and the expression simplifies to

$$\begin{equation} S_z \approx \frac{B}{2\sqrt{\mu }} \sqrt{\rho } \left(\frac{\rho _{\rm i}}{\rho _{\rm n}} v_{{\rm i},0} + v_{{\rm n},0}\right)^2. \end{equation} \tag{ 80 }$$

In the analysis of Vranjes et al. (2008) the initial disturbance is localized in the ionized fluid only. So Vranjes et al. (2008) took v_{n, 0} = 0. In this case Equation (80) becomes

$$\begin{equation} S_z \approx \frac{B}{2\sqrt{\mu }} \sqrt{\rho } \left(\frac{\rho _{\rm i}}{ \rho _{\rm n}} \right)^2v_{{\rm i},0}^2. \end{equation} \tag{ 81 }$$

Due to the presence of the factor (ρ_i/ρ_n)² the energy flux computed by Vranjes et al. (2008) in a weakly ionized plasma is much lower than the corresponding value in a fully ionized plasma (Equation (78)). The physical reason is that a significant portion of the wave energy is used to set into motion the neutral fluid, which is initially at rest in Vranjes et al. (2008).

On the contrary, Tsap et al. (2011) imposed $\hat{v}_{y} = v_{{\rm i},0}$, which can be satisfied only if v_{i, 0} = v_{n, 0} according to Equation (72). This means that Tsap et al. (2011) implicitly assumed that both ions and neutrals are initially perturbed with the same velocity, although this condition is not explicitly stated in the paper. The energy flux is then

$$\begin{equation} S_z \approx \frac{B}{2\sqrt{\mu }} \sqrt{\rho } \, v_{{\rm i},0}^2. \end{equation} \tag{ 82 }$$

The expression found by Tsap et al. (2011) is the same as Equation (78) for a fully ionized plasma. This is so because in Tsap et al. (2011) neutrals have initially the same velocity as ions and no wave energy needs to be used in setting neutrals into motion. Tsap et al. (2011) claimed that the expression of Vranjes et al. (2008) is incorrect. Here we clearly see that the discrepancy between Equations (81) and (82) is caused by a different choice of the initial conditions and that both expressions are indeed correct. The source of the discrepancy was expressing the energy flux in terms of the initial velocity of ions. This is a velocity definitely not easy to measure observationally. Instead, the velocity that an observer can measure is the velocity amplitude after the relaxation phase, $\hat{v}_{y}$. Then, Equation (77) should be used to estimate the energy flux.

Finally we recall that the present analysis, and so Equation (77), is valid in the case of Alfvén waves propagating in a uniform medium so that the energy flux is uniform. As commented in Section 5, the application of Equation (77) to the case of Alfvénic waves in the solar atmosphere must be done with caution since neither the density nor the wave velocity amplitude are uniform in the solar plasma (see Goossens et al. 2013).

In this paper we have studied theoretically Alfvén waves in a plasma composed of an ion–electron fluid and a neutral fluid interacting by means of neutral–ion collisions. This configuration is relevant in many astrophysical and laboratory plasmas. To keep our investigation as general as possible we have taken the neutral–ion collision frequency and the ionization degree as free parameters.

First, we have performed a normal mode analysis and have derived the dispersion relation of linear Alfvén waves. The dispersion relation agrees with the equations previously found by, e.g., Kulsrud & Pearce (1969), Pudritz (1990), Martin et al. (1997), Kumar & Roberts (2003), Zaqarashvili et al. (2011b), and Mouschovias et al. (2011). As in previous studies, we find that Alfvén waves are damped by neutral–ion collisions. The damping is most efficient when the wave frequency and the collision frequency are of the same order of magnitude. This conclusion applies to both standing and propagating waves. The effective Alfvén velocity of the plasma is modified by collisions, so that for high collision frequencies compared to the wave frequency the effective Alfvén velocity depends on the total density of the plasma and not on the ion density only (see also Kumar & Roberts 2003).

In addition, an important result is that when χ > 8, i.e., when the plasma is weakly ionized, there is a range of wavenumbers (or, equivalently, of wavelengths) for which oscillatory standing Alfvén waves are not possible. These cutoff wavenumbers are physical and are caused by a purely two-fluid effect (see Kulsrud & Pearce 1969; Pudritz 1990; Kamaya & Nishi 1998; Mouschovias et al. 2011). We investigated the physical reason for the existence of this cutoff region by analyzing the relative importance of the magnetic tension force and the neutral–ion friction force. We found that friction becomes the dominant force in the cutoff region, while tension is more important outside the cutoff interval. Hence, a disturbance in the magnetic field whose wavenumber is within the cutoff range decays due to collisions before the plasma is able to feel the restoring force of magnetic tension. As a consequence, oscillations of the magnetic field are effectively suppressed. The cutoff wavenumbers investigated here do not appear in the single-fluid approximation and are different from the unphysical cutoffs found in the single-fluid approximation (Zaqarashvili et al. 2012).

The solution of the initial-value problem is fully consistent with the normal modes analysis, and shows a growing strength of the interaction between ions and neutrals as the collision frequency increases. For high collision frequencies both ions and neutrals behave as a single fluid. Here a "high collision frequency" means that the neutral–ion collision frequency is at least an order of magnitude higher than the wave frequency. This condition is fulfilled in many astrophysical plasmas which makes the single-fluid approximation appropriate in those situations for the computations of periods/wavelengths and damping times/damping lengths. However, as explained above the single-fluid limit does not fully capture the details of the interaction between ions and neutrals and, for example, the existence of a range of cutoff wavenumbers is a two-fluid result.

As an example, we have considered Alfvén waves in a plasma with physical condition akin to those in the low solar atmosphere. As a matter of fact, due to the large values of the collision frequency this plasma could be studied using the single-fluid approximation instead of the more general two-fluid theory. In that respect, the effective Alfvén velocity is found to depend on the total density of the plasma. However, the presence of a certain range of cutoff wavelengths that can constrain the existence of oscillatory standing modes in the chromosphere is a result absent from those previous studies that use the limit of strong coupling (e.g., Haerendel 1992; De Pontieu & Haerendel 1998; Khodachenko et al. 2004; Leake et al. 2005). There are other astrophysical situations in which the ion–neutral coupling is not so strong as in the solar atmosphere and, therefore, a two-fluid treatment is needed. This may be the case of, e.g., protostellar discs (see a table with values of some physical parameters realistic of protostellar discs in Malyshkin & Zweibel 2011) and molecular clouds (see, e.g., Balsara 1996; Mouschovias et al. 2011).

In the present work we have restricted ourselves to Alfvén waves. Compressional magnetoacoustic waves have not been investigated. It is expected that both types of MHD waves are simultaneously present in a plasma. Magnetoacoustic waves in a two-fluid plasma will be investigated in a forthcoming work.

We acknowledge the support from MINECO and FEDER funds through grant AYA2011-22846 and from CAIB through the "grups competitius" scheme and FEDER funds. J.T. acknowledges support from MINECO through a Ramón y Cajal grant. R.S. thanks Ramon Oliver for reading the manuscript and for his constructive criticism.

The expressions of the coefficients of Equations (60) and (61) are as follows:

$$\begin{eqnarray} A_1 &=& \epsilon \frac{\left(\epsilon + \nu _{\mathrm{ni}}\right)\Gamma _{0, \rm i} + \chi \nu _{\mathrm{ni}}\Gamma _{0, \rm n}}{\omega _{\rm R}^2 + \left(\epsilon - \omega _{\rm I} \right)^2}, \end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray} A_2 &=& \frac{\left[ \omega _{\rm R}^2 + \omega _{\rm I}^2 - \epsilon \left(\nu _{\mathrm{ni}}+ 2 \omega _{\rm I} \right) \right]\Gamma _{0, \rm i} + \chi \nu _{\mathrm{ni}}\epsilon \Gamma _{0, \rm n}}{\omega _{\rm R}^2 + \left(\epsilon - \omega _{\rm I} \right)^2}, \end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray} A_3 &=& \frac{\omega _{\rm R}^2 \left(\epsilon + \nu _{\mathrm{ni}}+ \omega _{\rm I} \right) - \omega _{\rm I}\left(\epsilon - \omega _{\rm I} \right) \left(\nu _{\mathrm{ni}}+ \omega _{\rm I} \right)}{\omega _{\rm R}^2 + \left(\epsilon - \omega _{\rm I} \right)^2} \Gamma _{0, \rm i}\nonumber\\ && + \frac{\chi \nu _{\mathrm{ni}}\left(\omega _{\rm R}^2 + \omega _{\rm I}^2 - \epsilon \omega _{\rm I} \right)}{\omega _{\rm R} \left[ \omega _{\rm R}^2 + \left(\epsilon - \omega _{\rm I} \right)^2 \right]} \Gamma _{0, \rm n}, \end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray} C_1 &=& \frac{\nu _{\mathrm{ni}}\epsilon }{\epsilon + \nu _{\mathrm{ni}}} \frac{\left(\epsilon + \nu _{\mathrm{ni}}\right)\Gamma _{0, \rm i} + \chi \nu _{\mathrm{ni}}\Gamma _{0, \rm n}}{\omega _{\rm R}^2 + \left(\epsilon - \omega _{\rm I} \right)^2}, \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} C_2 &=& \left(1 + \frac{\chi \nu _{\mathrm{ni}}^3}{\left[ \omega _{\rm R}^2 + \left(\nu _{\mathrm{ni}}+ \omega _{\rm I} \right)^2 \right] \left(\epsilon + \nu _{\mathrm{ni}}\right)} \right) \Gamma _{0, \rm n}, \end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray} C_3 &=& {-} \frac{\epsilon \nu _{\mathrm{ni}}}{\omega _{\rm R}^2 + \left(\epsilon - \nu _{\mathrm{ni}}\right)^2} \Gamma _{0, \rm i} \nonumber\\ &&- \frac{\chi \nu _{\mathrm{ni}}^2 \left(\omega _{\rm R}^2 + \omega _{\rm I}^2 + \epsilon \nu _{\mathrm{ni}}\right)}{\left[ \omega _{\rm R}^2 + \left(\nu _{\mathrm{ni}}+ \omega _{\rm I} \right)^2 \right] \left[ \omega _{\rm R}^2 + \left(\epsilon - \omega _{\rm I} \right)^2\right]} \Gamma _{0, \rm n}, \end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray} C_4 &=& \frac{\nu _{\mathrm{ni}}\left(\omega _{\rm R}^2 + \omega _{\rm I}^2 - \epsilon \omega _{\rm I} \right)}{\omega _{\rm R}^2 + \left(\epsilon - \nu _{\mathrm{ni}}\right)^2} \Gamma _{0, \rm i} \nonumber\\ &&+ \frac{\chi \nu _{\mathrm{ni}}^2 \left[ \omega _{\rm R}^2 \left(\omega _{\rm I} - \epsilon +\nu _{\mathrm{ni}}\right) + \omega _{\rm I} \left(\nu _{\mathrm{ni}}+ \omega _{\rm I} \right) \left(\omega _{\rm I} - \epsilon \right) \right]}{\omega _{\rm R}\left[ \omega _{\rm R}^2 + \left(\nu _{\mathrm{ni}}+ \omega _{\rm I} \right)^2 \right] \left[ \omega _{\rm R}^2 + \left(\epsilon - \omega _{\rm I} \right)^2\right]} \Gamma _{0, \rm n}, \nonumber\\ \end{eqnarray} \tag{ A7 }$$

with the expressions of ω_R, ω_I, and given in Equations (21)–(23). In the expression of ω_R (Equation (21)) the + sign has to be taken.