© ESO 2020

1. Introduction

Ground- and space-based observations have confirmed the presence of oscillatory motions in structures of the solar atmosphere, such as coronal loops and prominences. They have been interpreted in terms of standing or propagating magnetohydrodynamic (MHD) waves, opening the door to the development of coronal and prominence seismology. This approach seeks to obtain information about physical conditions in these coronal structures by comparing observations and theoretical models of oscillations. Reviews on these topics are available in Roberts et al. (1984), Roberts (2000, 2008), Aschwanden (2004), Nakariakov & Verwichte (2005), Nakariakov et al. (2016), Oliver (1999, 2009), Erdélyi & Goossens (2011), Ballester (2015), Arregui et al. (2018). In the case of solar prominences, they display various types of oscillatory motions; Oliver (1999) classified these oscillations into two categories based on velocity amplitude: small-amplitude events with amplitudes in the range of 3 − 10 km s⁻¹, and large-amplitude events with amplitudes greater than 10 km s⁻¹. These two categories represent truly different phenomena: small-amplitude oscillations are, in general, not related to flare activity and only affect a small volume of the prominence. On the other hand, large-amplitude oscillations are often associated with an energetic event that sets the full prominence (or a large part of it) into an oscillatory state (Luna et al. 2018). Hence, it is likely more accurate to speak of local versus global oscillatory events in prominences.

Alfvén waves are of paramount importance in laboratory plasma physics and astrophysics. In the case of solar physics, these waves have been used to explain how the energy flows through and heats the solar atmosphere; however, direct detection of these waves has proved to be difficult. In recent years, claims regarding the detection of these elusive Alfvén waves in different layers and structures of the solar atmosphere, such as sunspot umbra, bright-points, and spicules, have been made (Tomczyk et al. 2007; Jess et al. 2009; Mathioudakis et al. 2013; Srivastava et al. 2017; Grant et al. 2018; Kohutova et al. 2020). However, most of these reports are based on indirect signatures that have been interpreted in terms of Alfvén waves. From a theoretical point of view, linear Alfvén waves have been thoroughly studied (Cramer 2001; Mathioudakis et al. 2013) under the assumption of fully ionized plasmas. However, in recent years, studies of Alfvén waves in partially ionized plasmas have been developed. The main reason for this increased study is that these plasmas constitute an essential ingredient of the solar atmosphere, since many layers, such as the chromosphere and photosphere, and structures, such as prominences and spicules, are made of partially ionized plasmas (Zaqarashvili et al. 2012; Soler et al. 2015a,b, 2016; Cally & Khomenko 2018, 2019; González-Morales et al. 2019; Khomenko & Cally 2019).

The damping of Alfvén waves has been studied by different authors using different approaches. For instance, Leake et al. (2005) studied the damping of linear Alfvén waves due to collisions between neutrals and ions in a Vernazza, Avrett, Loeser (VALC) model (Vernazza et al. 1981) for the chromosphere using analytical and numerical methods. The study allowed the authors to determine which frequencies are completely damped and which are unaffected and allow for an undamped propagation of the waves. Song (2011) studied the heat released by strongly damped Alfvén waves in the solar corona. In this case, a two-fluid model (plasma and neutrals) was used and the damping was produced by self-consistently considering the collisions among electrons, ions, and neutrals, as well as the interaction between charged particles and the electromagnetic field. The authors applied this model to the chromosphere and the results indicate that lower frequencies are nearly undamped and can propagate through the atmosphere, while higher frequencies are strongly damped at low altitudes. Threlfall et al. (2011) studied the damping of Alfvén waves in the ion-cyclotron frequencies when the Hall term is taken into account and when uniform and non-uniform plasmas are considered. In the case of uniform plasmas, they used an analytical approach, while for non-uniform plasmas a numerical approach was needed. Finally, Lazarian (2016) investigated the effect of turbulence on the damping of Alfvén waves, which may have different astrophysical applications ranging from the launching of stellar and galactic winds to the propagation of cosmic rays in galaxies and clusters of galaxies.

Nonlinear Alfvén waves in fully ionized plasmas have also been a subject of study. It is well known that the behavior of a linearly polarized Alfvén wave depends on wave amplitude, and that in the case of large amplitudes these waves generate density perturbations and motions along the magnetic field lines. In this scenario a self-interaction appears, which is due to the ponderomotive force. Hollweg (1971) solved the equations for a linearly polarized Alfvén wave propagating parallel to the direction of the magnetic field in a perfectly conducting fluid up to second order in the wave quantities. In the second order, longitudinal wave velocity and density fluctuations appear to be driven by gradients in the wave magnetic pressure. The main conclusions were: in the low β regime, large-amplitude linearly polarized Alfvén waves also have longitudinal motions and density fluctuations associated with them; and the Alfvén wave magnetic field and transverse velocity are not affected by second-order effects, but are affected by third-order terms, which were neglected. Rankin et al. (1994) and Tikhonchuk et al. (1995) investigated the nonlinear dynamics of standing shear Alfvén waves in cold magnetized plasmas. First, and in the limit β = 0, they obtained analytical approximations for second-order longitudinal velocity and density perturbations, corresponding to slow waves, as well as analytical approximations for third-order Alfvén wave quantities. In the second order, they found a secularly growing temporal behavior for longitudinal velocity and density perturbations; this produces a density enhancement at the velocity antinodes of Alfvén waves and a depletion near the nodes, with an angular frequency equal to 2k_zc_s, where c_s is the sound speed, and an effective wavenumber of 2k_z, where k_z is the longitudinal wavenumber of the Alfvén wave. When the assumption of β = 0 is removed, the secular effects on longitudinal velocity and density perturbations disappear and periodic behaviors are found. Finally, they derived third-order analytical approximations and used numerical solutions to the full set of nonlinear ideal MHD equations to test the validity of analytical results. Verwichte et al. (1999) studied the temporal evolution of Alfvénic pulses in cold homogeneous plasmas. The initial disturbance produces two pulses traveling in opposite directions; the ponderomotive force of the two pulses produces a static shock in longitudinal velocity at the initial position, while the traveling pulses form a shock described by a Cohen-Kulsrud equation. They found a good agreement between the derived analytical solutions and the numerical results. Terradas & Ofman (2004) studied the effect of large-amplitude waves in loops as a possible mechanism to produce density enhancements at loop tops due to the effect of the ponderomotive force of standing waves. They solved the nonlinear three-dimensional MHD equations in a flux tube configuration and found that for large initial perturbations a pressure imbalance appears along the loop, producing upflows from the footpoints. Therefore, the accumulation of mass at the loop top produces a significant density enhancement. Thurgood & McLaughlin (2013) numerically investigated the effects of ponderomotive force induced by Alfvén waves in inhomogeneous 2.5D MHD plasmas. They analysed the source terms in the nonlinear wave equations when the magnetic field and density are inhomogeneous; they conclude that the ponderomotive effects are induced by any Alfvén wave propagating in any medium and may affect the dynamics of energy transport and aspects of dissipation. Zheng et al. (2016) studied the effects of dissipative mechanisms, such as resistivity and viscosity, on nonlinear Alfvén wave trains. In weakly dissipative one-dimensional systems, the main effect is the damping of Alfvén waves and the heating of plasma; however, in the case of strong dissipation, the Alfvén wave train develops a damped soliton. In two-dimensional systems, phase mixing is present, which enhances damping and plasma heating. Finally, Lardner & Trehan (1991) studied the nonlinear evolution of Alfvén and magnetoacoustic waves in a low-density plasma with a strong magnetic field, using a very different approach based on the Chew-Goldberger-Low approximation. In the reported research on nonlinear Alfvén waves, the plasma has always been considered as fully ionized and treated as a single fluid. This assumption of fully ionized plasma is valid, for instance, for the solar corona; however, as we have previously pointed out, the plasma is only partially ionized in some layers and structures of the solar atmosphere. On the other hand, except in one case (Zheng et al. 2016), dissipative effects have not been taken into account, and ideal or resistive MHD equations have been used. Therefore, it is of great interest to study the behavior of nonlinear Alfvén waves in partially ionized plasmas; this approach was taken by Martínez-Gómez et al. (2018), who, using a nonlinear multi-fluid code, investigated the behavior of high-frequency nonlinear waves in a partially ionized plasma, paying special attention to the heating caused by ion-neutral collisions.

The aim of this paper is to study the temporal behavior of nonlinear Alfvén waves in a partially ionized plasma, whose physical conditions are those of quiescent solar prominences, when ambipolar diffusion is taken into account. For this study, we consider our plasma as a single fluid and, apart from ambipolar diffusion, we have also included radiative losses and thermal conduction as dissipative mechanisms for the field-aligned motions and density perturbations excited in the plasma by the ponderomotive force. Furthermore, third-order Alfvén waves, which are influenced by these second-order perturbations which in turn are influenced by thermal effects, have also been studied.

The layout of the paper is as follows: First we introduce single-fluid MHD equations, an equilibrium model and dissipative mechanisms; next, we derive equations for parallel propagating waves and apply the regular perturbations method, obtaining first-, second-, and third-order equations which have been solved for different cases; finally, we obtain numerical solutions for the full set of nonlinear MHD equations and compare the results with those from the regular perturbations approach.

2. Basic equations

2.1. Governing equations in the single-fluid approximation

We study the dynamics of partially ionized plasmas within the framework of the single-fluid approximation. In short, the single-fluid approximation assumes a strong coupling between ions, electrons, and neutrals so that all the species effectively behave as one fluid. In this approximation, the basic MHD equations are written in terms of total quantities, while the effect of the interactions between the various species remains in the form of several nonideal terms. For instance, the influence of ion-neutral collisions is present through the so-called ambipolar diffusion effect, which acts as a diffusive mechanism for the magnetic field. The single-fluid approximation is appropriate when studying MHD waves in partially ionized plasmas as long as the wave frequency remains lower than the ion-neutral collision frequency. The values of the collision frequencies in the solar chromosphere can be checked, for example, in Fig. 11 of Ballester et al. (2018); in solar prominences, the collision frequencies are similar to those found in the upper chromosphere. In the following expressions, MKS units have been used.

The basic single-fluid equations used in this work are:

$$\begin{array}{cc}& \frac{\mathrm{D}\rho }{\mathrm{D}t}=-\rho \mathrm{\nabla }·\mathit{v},\hfill \end{array}$$(1)

$$\begin{array}{cc}& \rho \frac{\mathrm{D}\mathit{v}}{\mathrm{D}t}=-\mathrm{\nabla }p+\frac{1}{\mu }(\mathrm{\nabla }\times \mathit{B})\times \mathit{B},\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial \mathit{B}}{\partial t}=\mathrm{\nabla }\times (\mathit{v}\times \mathit{B})-\mathrm{\nabla }\times (\eta \mathrm{\nabla }\times \mathit{B})\hfill \\ \hfill & +\mathrm{\nabla }\times \{{\eta }_{\mathrm{A}}[(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}]\times \mathit{B}\},\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\mathrm{D}p}{\mathrm{D}t}=-\gamma p\mathrm{\nabla }·\mathit{v}+(\gamma -1)\mathcal{L},\hfill \end{array}$$(4)

$$\begin{array}{cc}& p=\frac{\rho RT}{\stackrel{\sim }{\mu }},\hfill \end{array}$$(5)

where $\frac{\text{D}}{\text{D}t}=\frac{\partial }{\partial t}+\mathit{v}·\mathrm{\nabla }$ denotes the material or total derivative, ρ is the mass density, p is the thermal pressure, T is the temperature, v is the velocity vector, B is the magnetic field vector, γ is the adiabatic index, μ the magnetic permeability, η and η_A are the coefficients of Ohmic and ambipolar diffusion, respectively, ℒ represents the net effect of all the sources and sinks of energy, R is the gas constant, and $\stackrel{\sim }{\mu }$ is the mean molecular weight. Equations (1)–(5) are the continuity, momentum, induction, energy, and state equations, respectively.

2.2. Ionization state

The ionization state of the plasma is determined by the balance between ionization and recombination. The general approach to computing the ionization state is to take into account all the ionization and recombination processes and to determine the time-dependent values of the ionized and neutral fractions. Here, to estimate the ionization degree, we follow Heinzel et al. (2015), who used one-dimensional non-LTE radiative transfer models (Heinzel et al. 2014) to determine the ionization degree in several prominence slabs. In particular, these authors provide tables for the ionization degree for different temperatures and pressures at the prominence. Heinzel et al. (2015) consider a prominence plasma composed of hydrogen and helium, whose abundance is 10% and which is fully neutral. Then, the ionization degree i is defined as $i=\frac{n_{\mathrm{e}}}{n_{\mathrm{H}}}$, where n_e is the electron density number and n_H is the total hydrogen density number (n_H I + n_p), using the following subscripts: e for electrons, p for protons, H for neutral hydrogen plus protons, HI for neutral hydrogen, and HeI for neutral helium; the total particle density number is N = n_H + n_He I +n_e. Then, taking into account that n_He| I = 0.1n_H, N can be written as

$$\begin{array}{c}\hfill N=n_{\mathrm{e}}(1+\frac{1.1}{i})\end{array}$$(6)

and the gas pressure, p = Nk_BT, is given as

$$\begin{array}{c}\hfill p=n_{\mathrm{e}}(1+\frac{1.1}{i})k_{\mathrm{B}}T\end{array}$$(7)

or, using the mass density, as

$$\begin{array}{c}\hfill p=\frac{\rho k_{\mathrm{B}}T(i+1.1)}{1.4H},\end{array}$$(8)

where k_B is the Boltzmann constant and H is the atomic mass unit. Using Table 1 in Heinzel et al. (2015) and the Fit function from the Mathematica symbolic package, we performed a polynomial fit, up to the third order in pressure p and temperature T as well as product terms and interactions, of the ionization degree:

$$\begin{array}{c}\hfill i=a+bp+c{p}^2+d{p}^3+eT+f{T}^2+g{T}^3+hpT+j{p}^{2}T+k{T}^{2}p,\end{array}$$(9)

where a, b, c, d, e, f, g, h, j, and k are the coefficients of the fitted function (see Appendix A). Then, using Eq. (9) together with the temperature and pressure values in Heinzel et al. (2015), we can compute the fitted ionization degree. The goodness of our surface fitting is assessed by computing the determination coefficient R²:

$$\begin{array}{c}\hfill R^2=1-\frac{{\displaystyle \sum _i{(i_i-{\widehat{i}}_i)}^2}}{{\displaystyle \sum _i{(i_i-{\overline{i}}_i)}^2}},\end{array}$$(10)

where i_i corresponds to the values of the ionization degree in Heinzel et al. (2015), ${\widehat{i}}_i$ corresponds to the fitted ionization degree using Eq. (9), and ${\overline{i}}_i$ corresponds to the mean of the values of the ionization degree in Heinzel et al. (2015). Using Eq. (10), the obtained value for R² is 0.98 or 98%. Figure 1 compares the fitted surface, described by Eq. (9), to data in Table 1 of Heinzel et al. (2015), visually showing the goodness of the fit.

Fig. 1. Ionization degree i versus temperature T and pressure p. The blue dots correspond to data in Table 1 in Heinzel et al. (2015); the surface represents the fitting described by Eq. (9).

Next, in order to determine the total particle density number we need to know the plasma pressure; from Eq. (8), the ionization degree can be written as

$$\begin{array}{c}\hfill i=\frac{p1.4H}{\rho k_{\mathrm{B}}T}-1.1.\end{array}$$(11)

Then, once a value for the density was assumed, and using the different temperature values in the Heinzel et al. (2015) table, we solved Eqs. (9) and (11) together. When the pressure is determined, the total particle density number is obtained from

$$\begin{array}{c}\hfill N=\frac{p}{k_{\mathrm{B}}T}\end{array}$$(12)

and the rest of particle density numbers can be determined from

$$\begin{array}{cc}& n_{\mathrm{e}}=n_{\mathrm{p}}=\frac{N}{1+\frac{1.1}{i}},\hfill \end{array}$$(13)

$$\begin{array}{cc}& n_{\mathrm{H}}=\frac{n_{\mathrm{e}}}{i},\hfill \end{array}$$(14)

$$\begin{array}{cc}& n_{\text{He}\text{I}}=N-n_{\mathrm{H}}-n_{\mathrm{e}},\hfill \end{array}$$(15)

$$\begin{array}{cc}& n_{\text{H}\text{I}}=n_{\mathrm{H}}-n_{\mathrm{e}}.\hfill \end{array}$$(16)

Finally, the mean molecular weight can be computed as

$$\begin{array}{c}\hfill \stackrel{\sim }{\mu }=\frac{\overline{m}}{H}=\frac{\frac{{\sum }_{\beta }{n}_{\beta }{m}_{\beta }}{{\sum }_{\beta }{n}_{\beta }}}{H}=\frac{(n_{\mathrm{p}}{m}_{\mathrm{p}}+n_{\mathrm{HI}}{m}_{\mathrm{HI}}+n_{\mathrm{HeI}}{m}_{\mathrm{HeI}}+n_{\mathrm{e}}{m}_{\mathrm{e}})}{(n_{\mathrm{p}}+n_{\mathrm{HI}}+n_{\mathrm{HeI}}+n_{\mathrm{e}})H},\end{array}$$(17)

where $\overline{m}$ is the mean mass per particle.

2.3. Dissipative mechanisms

The induction equation (Eq. (3)) contains two magnetic diffusion terms in addition to the ideal convective term. The second term on the right-hand side is Ohmic diffusion, which is caused by electron collisions. The coefficient of Ohmic diffusion is given by

$$\begin{array}{c}\hfill \eta =\frac{{\alpha }_{\mathrm{e}}}{{\mu }_0{e}^2{n}_{\mathrm{e}}^2},\end{array}$$(18)

where e is the electron charge. The total friction coefficient of species β with other species is

$$\begin{array}{c}\hfill {\alpha }_{\beta }={\displaystyle \sum _{\beta \prime \ne \beta }{\alpha }_{\beta \beta \prime }.}\end{array}$$(19)

Therefore, in our case, the total friction coefficient of electrons α_e with protons, neutral hydrogen, and neutral helium is

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{e}}={\alpha }_{\mathrm{ep}}+{\alpha }_{e\text{H}\text{I}}+{\alpha }_{e\text{He}\text{I}}.\end{array}$$(20)

The friction coefficient for collisions between two charged species q and q′ (Spitzer 1962; Braginskii 1965) is

$$\begin{array}{c}\hfill {\alpha }_{qq\prime }=\frac{n_{\mathrm{q}}{n}_{\mathrm{q}\prime }{e}^{4}ln{\mathrm{\Lambda }}_{qq\prime }}{6\pi \sqrt{2\pi }{ϵ}_0^2{m}_{\mathrm{qq}\prime }{(k_{\mathrm{B}}T/m_{\mathrm{qq}\prime })}^{3/2}},\end{array}$$(21)

where m_qq′ = m_qm_q′/(m_q+m_q′) is the reduced mass, ϵ₀ is the permittivity of free space, and lnΛ_qq′ is the Coulomb’s logarithm (Spitzer 1962; Vranjes & Krstic 2013). This logarithm is given by

$$\begin{array}{c}\hfill ln{\mathrm{\Lambda }}_{qq\prime }=ln\left(\frac{24\pi {ϵ}_0^{1.5}{k}^{1.5}{T}^{1.5}}{e^3\sqrt{n_q+n_{q\prime }}}\right),\end{array}$$(22)

while the friction coefficient for collisions between a charged or neutral species β = e, p, H I, He I, and neutral species n = H I, He I is

$$\begin{array}{c}\hfill {\alpha }_{\beta n}=m_{\beta }{n}_n{m}_{\beta n}{\left[\frac{8kT}{\pi m_{\beta n}}\right]}^{1/2}{\sigma }_{\beta n},\end{array}$$(23)

where σ_βn is the collision cross section (Vranjes & Krstic 2013).

The third term on the right-hand side of Eq. (3) is ambipolar diffusion. This term contains the effect of ion-neutral collisions in the single-fluid approximation. The coefficient of ambipolar diffusion in a hydrogen-helium plasma (Zaqarashvili et al. 2013) is

$$\begin{array}{c}\hfill {\eta }_{\mathrm{A}}=\frac{{\xi }_{\text{H}\text{I}}^2{\alpha }_{\text{He}\text{I}}+{\xi }_{\text{He}\text{I}}^2{\alpha }_{\text{H}\text{I}}+2{\xi }_{\text{H}\text{I}}{\xi }_{\text{He}\text{I}}{\alpha }_{\text{H}\text{I}\text{He}\text{I}}}{{\mu }_0({\alpha }_{\text{H}\text{I}}{\alpha }_{\text{He}\text{I}}-{\alpha }_{\text{H}\text{I}\text{He}\text{I}}^2)},\end{array}$$(24)

with the friction coefficients for α_HeI and α_HI given by

$$\begin{array}{cc}& {\alpha }_{\text{He}\text{I}}={\alpha }_{\text{He}\text{I}p}+{\alpha }_{\text{He}\text{I}\text{H}\text{I}}+{\alpha }_{\text{He}\text{I}e}\hfill \end{array}$$(25)

$$\begin{array}{cc}& {\alpha }_{\text{H}\text{I}}={\alpha }_{\text{H}\text{I}p}+{\alpha }_{\text{H}\text{I}\text{He}\text{I}}+{\alpha }_{\text{H}\text{I}e}.\hfill \end{array}$$(26)

The effect of all the sources and sinks of energy is included in the function ℒ on the right-hand side of the energy equation (Eq. (4)). The expression of ℒ is

$$\begin{array}{c}\hfill \mathcal{L}=\mathrm{\nabla }·\left(\kappa \mathrm{\nabla }T\right)-\rho L(\rho ,T)+\mathbf{j}·{\mathbf{E}}^{\ast }+\rho h,\end{array}$$(27)

where the various terms on the right-hand side are the heat flux due to thermal conduction ∇ ⋅ (κ∇T). Here, κ is the thermal conductivity tensor and L(ρ,T) is the radiative loss function. The j ⋅ E^* is the generalized Joule heating, where $\mathbf{j}$ is the current density and E^* is the effective electric field. The h represents an additional and unspecified external constant heating per unit volume. The thermal conduction in a magnetized plasma is strongly anisotropic. For convenience, we split the thermal conduction term in Eq. (27) into its components that are parallel and perpendicular to the magnetic field direction, namely

$$\begin{array}{c}\hfill \mathrm{\nabla }·\left(\kappa \mathrm{\nabla }T\right)=\mathit{B}·\mathrm{\nabla }\left({\kappa }_{\Vert }\frac{\mathit{B}·\mathrm{\nabla }T}{{\left|\mathit{B}\right|}^2}\right)+\mathrm{\nabla }·\left({\kappa }_{\perp }\frac{\mathit{B}\times (\mathrm{\nabla }T\times \mathit{B})}{{\left|\mathit{B}\right|}^2}\right),\end{array}$$(28)

where κ_∥ and κ_⊥ are the parallel and perpendicular conductivities, respectively, with respect to the magnetic field direction, and |B|² is the square of the modulus of the magnetic field vector. In a fully ionized plasma, κ_∥ is essentially determined by the electron conductivity, while κ_⊥ is mainly determined by ions. In a partially ionized plasma, the isotropic conductivity of neutrals has to be added to both parallel and perpendicular components. Thus, the parallel and perpendicular components of thermal conductivity are approximated by κ_∥ ≈ κ_e, ∥ + κ_n and κ_⊥ ≈ κ_n, where κ_e, ∥ and κ_n are the parallel electron thermal conductivity and the total neutral thermal conductivity, respectively, as given by (Braginskii 1965; Meier 2011; Soler et al. 2015a):

$$\begin{array}{cc}& {\kappa }_{\mathrm{e},\Vert }=3.2\frac{n_{\mathrm{e}}^2{k}_{\mathrm{B}}^{2}T}{{\alpha }_{\mathrm{e},\mathrm{tot}}}\hfill \end{array}$$(29)

$$\begin{array}{cc}& {\kappa }_{\mathrm{n}}=\frac{5}{3}(\frac{n_{\mathrm{HI}}^2}{{\alpha }_{\mathrm{HI},\mathrm{tot}}}+\frac{n_{\mathrm{HeI}}^2}{{\alpha }_{\mathrm{HeI},\mathrm{tot}}})k_{\mathrm{B}}^{2}T,\hfill \end{array}$$(30)

where α_β, tot is the total friction coefficient of species β, given by α_β, tot = α_β + α_ββ; this includes the term α_ββ, which accounts for self-collisions.

The radiative loss function, L(ρ,T), accounts for radiative losses due to plasma cooling. Determining the plasma radiative losses as a function of temperature and density is a difficult task that requires complicated numerical solutions of the radiative transfer problem. The radiative loss rate depends, for example, on the completeness of the atomic model used for the calculation, on the atomic processes included, on the ionization equilibrium, and on the assumed element abundance, among other factors. The solution of the radiative transfer problem is beyond the scope of the present paper. A frequent alternative approach to account for the plasma radiative losses is to use semi-empirical parametrizations of the radiative loss function that are widely available in the literature (Rosner et al. 1978; Milne et al. 1979; Dahlburg & Mariska 1988; Klimchuk & Cargill 2001; Carbonell et al. 2004; Landi et al. 2012). These parametrizations enable us to incorporate radiative losses in a consistent way without needing to solve the full radiative transfer problem. The expression of the radiative loss function we used is

$$\begin{array}{c}\hfill L(\rho ,T)=\rho {\chi }^{\ast }{T}^{\alpha },\end{array}$$(31)

where χ^* and α are piecewise constants dependant on the temperature, and we use the parametrization of χ^* and α given in Hildner (1974) (see Table 1). The inconvenience of this approach is that the semi-empirical radiative loss function is obtained under the assumption of optically thin plasma, while the cool partially ionized plasmas of interest here may not completely satisfy this condition. Owing to finite optical thickness, the actual radiative losses of the plasma would probably be reduced to some degree compared to the optically thin calculation.

Finally, the generalized Joule heating j ⋅ E^* takes into account plasma heating due to the dissipation of electric currents by both Ohmic and ambipolar diffusion. The expression of j ⋅ E^* is

$$\begin{array}{c}\hfill \mathbf{j}·{\mathbf{E}}^{\ast }=\mu \eta {\left|{\mathbf{j}}_{\Vert }\right|}^2+\mu {\eta }_{\mathrm{C}}{\left|{\mathbf{j}}_{\perp }\right|}^2,\end{array}$$(32)

where η_C is the so-called Cowling’s (or total) diffusivity given by

$$\begin{array}{c}\hfill {\eta }_{\mathrm{C}}=\eta +{\left|\mathit{B}\right|}^2{\eta }_{\mathrm{A}}\end{array}$$(33)

and j_∥ and j_⊥ are the components of the current density parallel and perpendicular to the magnetic field, respectively, that can be computed as

$$\begin{array}{cc}& {\mathbf{j}}_{\Vert }=\frac{1}{\mu }\frac{[(\mathrm{\nabla }\times \mathit{B})·\mathit{B}]\mathit{B}}{{\left|\mathit{B}\right|}^2}\hfill \end{array}$$(34)

$$\begin{array}{cc}\hfill & {\mathbf{j}}_{\perp }=\frac{1}{\mu }\frac{\mathit{B}\times [(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}]}{{\left|\mathit{B}\right|}^2}·\hfill \end{array}$$(35)

Thus, Ohmic diffusion is responsible for the dissipation of parallel currents, while Cowling’s diffusion, that is to say the joint effect of Ohmic and ambipolar diffusion, is responsible for the dissipation of perpendicular currents.

2.4. Equations for parallel-propagating waves

We use a Cartesian coordinate system and assume that the plasma properties vary on the z-direction only, whereas x and y are ignorable directions. Then, we can conveniently rotate the coordinate system for v and B to lie on the zy-plane with no loss of generality. Hence, v = (0, v_y, v_z) and B = (0, B_y, B_z). The solenoidal condition ∇ ⋅ B = 0 imposes that ∂B_z/∂z = 0, while from the z-component of Eq. (3) we get ∂B_z/∂t = 0. Therefore, B_z is a constant in both space and time. In this one-dimensional case, Eqs. (1)–(4) become

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}=-\frac{\partial }{\partial z}\left(\rho v_z\right),\hfill \end{array}$$(36)

$$\begin{array}{cc}& \frac{\partial v_z}{\partial t}=-v_z\frac{\partial v_z}{\partial z}-\frac{1}{\rho }\frac{\partial p}{\partial z}-\frac{B_y}{\mu \rho }\frac{\partial B_y}{\partial z},\hfill \end{array}$$(37)

$$\begin{array}{cc}& \frac{\partial v_y}{\partial t}=-v_z\frac{\partial v_y}{\partial z}+\frac{B_z}{\mu \rho }\frac{\partial B_y}{\partial z},\hfill \end{array}$$(38)

$$\begin{array}{cc}& \frac{\partial B_y}{\partial t}=B_z\frac{\partial v_y}{\partial z}-\frac{\partial }{\partial z}\left(B_y{v}_z\right)+\frac{\partial }{\partial z}\left({\eta }_{\mathrm{C}}\frac{\partial B_y}{\partial z}\right),\hfill \end{array}$$(39)

$$\begin{array}{cc}& \frac{\partial p}{\partial t}=-v_z\frac{\partial p}{\partial z}-\gamma p\frac{\partial v_z}{\partial z}\hfill \\ \hfill & +(\gamma -1)[\frac{1}{\mu }{\eta }_{\mathrm{C}}{\left(\frac{\partial B_y}{\partial z}\right)}^2+\frac{\partial }{\partial z}\left(\frac{B_z^2{\kappa }_{\Vert }+B_y^2{\kappa }_{\perp }}{B_z^2+B_y^2}\frac{\partial T}{\partial x}\right)]\hfill \\ \hfill & -(\gamma -1)({\rho }^2{\chi }^{\ast }{T}^{\alpha }-\rho h),\hfill \end{array}$$(40)

3. Regular perturbations approach

Before exploring the numerical solutions to Eqs. (36)–(40), we perform an approximate analytic study of the solutions in the limit of small amplitudes in this section. To do so, we consider small-amplitude Alfvén waves superimposed on a homogeneous and constant equilibrium state. The purpose of this section is to establish an analytic basis that can be compared to the fully numerical simulations performed later. We define the parameter ϵ ≡ v₀/c_A, 0, where v₀ is the Alfvén wave initial velocity amplitude and c_A, 0 is the equilibrium Alfvén speed. Then, assuming that ϵ ≪ 1, we write

$$\begin{array}{cc}& \rho ={\rho }_0+{ϵ}^2\rho \prime ,\hfill \end{array}$$(41)

$$\begin{array}{cc}& p=p_0+{ϵ}^{2}p\prime ,\hfill \end{array}$$(42)

$$\begin{array}{cc}& T=T_0+{ϵ}^{2}T\prime ,\hfill \end{array}$$(43)

$$\begin{array}{cc}& v_z={ϵ}^2{v}_z\prime ,\hfill \end{array}$$(44)

$$\begin{array}{cc}& v_y=ϵv_y\prime +{ϵ}^{3}v{\prime }_{3y},\hfill \end{array}$$(45)

$$\begin{array}{cc}& B_z=B_0,\hfill \end{array}$$(46)

$$\begin{array}{cc}& B_y=ϵB_y\prime +{ϵ}^{3}B{\prime }_{3y},\hfill \end{array}$$(47)

where the subscript 0 denotes a homogeneous and constant equilibrium quantity and the prime ′ denotes a perturbation. Since, initially, we are interested in Alfvén waves, the perturbations of the perpendicular components of velocity and magnetic field are assumed to be of first order in ϵ with a third-order correction, while the perturbations of the remaining quantities are assumed to be second order in ϵ. Contributions with higher orders in ϵ are neglected. We substitute these quantities in Eqs. (36)–(40) and separate the various terms according to their order with respect to ϵ. From the zeroth order equations in ϵ, we find that the unspecified external heating source must be $h={\rho }_0{\chi }^{*}{T}_0^{\alpha }$ to compensate for radiative losses in order for the equilibrium to be satisfied. The first-order equations in ϵ govern the behavior of $v_y^{\prime }$ and $B_y^{\prime }$, and so they describe linear Alfvén waves. The second-order equations in ϵ involve $v_z^{\prime }$, p′, ρ′, and T′, so that they describe the generation of slow magnetoacoustic waves due to the nonlinear coupling with the Alfvén waves. The equations for higher orders in ϵ represent nonlinear corrections on the Alfvén and slow waves. In this approximate study, we restrict ourselves to the first-, second-, and third-order equations in ϵ.

3.1. First-order equations

In order to write all the equations under study in dimensionless forms, we introduce the following dimensionless quantities:

$$\begin{array}{cc}& \overline{v}\prime =\frac{v\prime }{c_{\mathrm{A},0}},{\overline{k}}_z=k_{z}L,\overline{z}=\frac{z}{L},\overline{t}=\frac{t}{\tau },{\overline{\eta }}_{C,0}=\frac{{\eta }_{C,0}}{c_{\mathrm{A},0}L},\hfill \\ \hfill & \overline{B}\prime =\frac{B\prime }{B_0},\overline{\rho }\prime =\frac{\rho \prime }{{\rho }_0},{\overline{c}}_{A,0}=\frac{c_{\mathrm{A},0}}{c_{\mathrm{A},0}}=1,\hfill \end{array}$$(48)

where L corresponds to half the size of the spatial domain under consideration and τ is a characteristic timescale, which are related through the equilibrium Alfvén speed c_A,0 = L/τ. Hereafter we drop bars for the sake of simplicity.

The first-order equations in ϵ are:

$$\begin{array}{cc}& \frac{\partial v_y\prime }{\partial t}=\frac{\partial B_y\prime }{\partial z}\hfill \end{array}$$(49)

$$\begin{array}{cc}& \frac{\partial B_y\prime }{\partial t}=\frac{\partial v_y\prime }{\partial z}+{\eta }_{\mathrm{C},0}\frac{{\partial }^2{B}_y\prime }{\partial z^2}·\hfill \end{array}$$(50)

These two equations involve the components of v and B transverse to the equilibrium magnetic field and can appropriately be combined into an equation for $v_y^{\prime }$ only, namely

$$\begin{array}{c}\hfill \frac{{\partial }^2{v}_y\prime }{\partial t^2}-\frac{{\partial }^2{v}_y\prime }{\partial z^2}-{\eta }_{\mathrm{C},0}\frac{{\partial }^3{v}_y\prime }{\partial z^2\partial t}=0,\end{array}$$(51)

where η_C,0 is Cowling’s diffusivity computed at the equilibrium state. Equation (51) governs linear Alfvén waves that are damped by ambipolar and Ohm’s diffusion, and it can also be shown that $B_y^{\prime }$ satisfies a similar equation. Since we are interested in an analytical solution to Eq. (51) describing damped linear standing Alfvén waves, we use the separation of variables method. We assume that $v_y^{\prime }$(z, t) can be written as

$$\begin{array}{c}\hfill v_y\prime (z,t)=u(z)g(t).\end{array}$$(52)

Then, substituting this expression in Eq. (51) and using $k_z^2$ as the separation constant, we obtain

$$\begin{array}{cc}& \frac{{\mathrm{d}}^{2}g(t)}{\mathrm{d}{t}^2}+k_z^2{\eta }_{\mathrm{C},0}\frac{\mathrm{d}g(t)}{\mathrm{d}t}+k_z^{2}g(t)=0\hfill \end{array}$$(53)

$$\begin{array}{cc}& \frac{{\mathrm{d}}^{2}u(z)}{\mathrm{d}{z}^2}+k_z^{2}u(z)=0,\hfill \end{array}$$(54)

where the separation constant k_z plays the role of the wavenumber along the magnetic field direction. Assuming the stationary state of wave propagation, the solution to Eq. (53) can be obtained assuming that g(t) behaves harmonically as e^iωt, which provides a dispersion relation of

$$\begin{array}{c}\hfill {\omega }^2-i{k}_z^2{\eta }_{\mathrm{C},0}\omega -k_z^2=0,\end{array}$$(55)

whose solution for ω is

$$\begin{array}{c}\hfill \omega =\frac{i{\eta }_{\mathrm{C},0}{k}_z^2}{2}\pm \frac{\sqrt{4k_z^2-{\eta }_{\mathrm{C},0}^2{k}_z^4}}{2},\end{array}$$(56)

a complex angular frequency, ω = ω_r + iω_i, with

$$\begin{array}{cc}& {\omega }_r=\pm \frac{\sqrt{4k_z^2-{\eta }_{\mathrm{C},0}^2{k}_z^4}}{2}\hfill \end{array}$$(57)

$$\begin{array}{cc}\hfill & {\omega }_i=\frac{{\eta }_{\mathrm{C},0}{k}_z^2}{2}·\hfill \end{array}$$(58)

Then, g(t) can be written as

$$\begin{array}{c}\hfill g(t)=e^{-\frac{1}{2}{k}_z^2{\eta }_{\mathrm{C},0}t}(C_{1}cos{\omega }_{r}t+C_{2}sin{\omega }_{r}t).\end{array}$$(59)

On the other hand, the solution to Eq. (54) is

$$\begin{array}{c}\hfill u(z)=C_{3}cos{k}_{z}z+C_{4}sin{k}_{z}z.\end{array}$$(60)

Then, $v_y^{\prime }$(z, t) can be written as

$$\begin{array}{cc}& v_y\prime (z,t)=e^{-\frac{1}{2}{k}_z^2{\eta }_{\mathrm{C},0}t}(C_{1}cos{\omega }_{r}t+C_{2}sin{\omega }_{r}t)\hfill \\ \hfill & \times (C_{3}cos{k}_{z}z+C_{4}sin{k}_{z}z)\hfill \end{array}$$(61)

Now, in order to determine the four constants of integration, we impose the following conditions representative of a standing oscillation with wavelength equal to 4L:

$$\begin{array}{cc}& v_y\prime (z=-1,t)=0,\hfill \end{array}$$(62)

$$\begin{array}{cc}& v_y\prime (z=1,t)=0,\hfill \end{array}$$(63)

$$\begin{array}{cc}& v_y\prime (z,0)=v_{0}cos{k}_{z}z,\hfill \end{array}$$(64)

$$\begin{array}{cc}& {\left(\frac{\partial v_y\prime (z,t)}{\partial t}\right)}_{t=0}=0,\hfill \end{array}$$(65)

where the longitudinal wavenumber $k_z=\frac{\pi }{2}$, and v₀ is the initial amplitude. Then, using the above conditions, the final solution for the transverse velocity perturbation is

$$\begin{array}{cc}& v_y\prime (z,t)=v_{0}exp[-\frac{1}{2}{k}_z^2{\eta }_{\mathrm{C},0}t](cos{\omega }_{r}t+\frac{k_z^2{\eta }_{\mathrm{C},0}}{2{\omega }_r}sin{\omega }_{r}t)\times cos{k}_{z}z\hfill \end{array}$$(66)

with the period and damping time given by

$$\begin{array}{cc}& P=\frac{2\pi }{{\omega }_r}=\frac{4\pi }{\sqrt{4k_z^2-{\eta }_{\mathrm{C},0}^2{k}_z^4}}\hfill \end{array}$$(67)

$$\begin{array}{cc}& {\tau }_D=\frac{1}{{\omega }_i}=\frac{2}{k_z^2{\eta }_{\mathrm{C},0}}·\hfill \end{array}$$(68)

Once the expression for the perturbed velocity amplitude is known, the solution for the perturbed magnetic field $B_y^{\prime }$ can be obtained by integrating

$$\begin{array}{c}\hfill \frac{\partial v_y\prime }{\partial t}=\frac{\partial B_y\prime }{\partial z}\end{array}$$(69)

to obtain

$$\begin{array}{c}\hfill B_y\prime (z,t)=\frac{v_{0}exp[-\frac{1}{2}{k}_z^2{\eta }_{C,0}t](k_z^4{\eta }_{\mathrm{C},0}^2+4{\omega }_r^2)sin{k}_{z}zsin{\omega }_{r}t}{4k_z{\omega }_r}.\end{array}$$(70)

In order to plot the temporal behavior of first-order perturbed velocity and magnetic fields, we impose the dimensionless initial amplitude of the velocity perturbation v₀ = 0.15 and compute other dimensionless quantities such as η_C,0 and ω_r. To this end, we assume values for the magnetic field, density, and temperature typical of quiescent prominences B₀ = 5 × 10⁻⁴ T, ρ₀ = 5 × 10⁻¹¹ kg ⋅ m⁻³, and T = 12 000 K. Then, using these values, the Alfvén speed, c_A,0, is 63.08 km s⁻¹ and the sound speed, c_s,0, is 15.08 km s⁻¹. Furthermore, in order to remain in the low-frequency regime and to satisfy the requirements of single-fluid approximation, we take L = 10⁶ m as our characteristic length and τ = 15.85 s is our characteristic timescale; these values give a period for the Alfvén waves of 64.35 s and a frequency of 0.015 Hz. Finally, using Eq. (48), our dimensionless quantities are η_C,0 = 0.0023 and ω_r = 1.57. The temporal behavior of the linear Alfvén velocity and magnetic field perturbations can be seen in Fig. 2, which shows the damping of the linear Alfvén wave produced by ambipolar diffusion.

Fig. 2. Dimensionless linear Alfvén wave velocity (red) and magnetic field perturbations (blue) versus dimensionless time. The velocity amplitude has been computed at z = 0, while the magnetic field amplitude has been computed at z = 1 (T = 12 000, v0 = 0.15, kz = π/2, ηC,0 = 0.0023).

3.2. Second-order equations

The second-order equation in ϵ corresponding to the energy equation (Eq. 40) is

$$\begin{array}{cc}\hfill \frac{\partial p\prime }{\partial t}& =-\gamma p_0\frac{\partial v_z\prime }{\partial z}+(\gamma -1)[{\kappa }_{\Vert ,0}\frac{{\partial }^{2}T\prime }{\partial z^2}-{\rho }_0(L_{\rho }\rho \prime +L_{T}T\prime )]\hfill \\ \hfill & +(\gamma -1)\left[\frac{1}{{\mu }_0}{\eta }_{\mathrm{C},0}{\left(\frac{\partial B_y\prime }{\partial z}\right)}^2\right],\hfill \end{array}$$(71)

where κ_∥, 0 denotes the value of κ_∥ computed with the equilibrium quantities (see Sect. 2.3), and L_ρ and L_T are the partial derivatives of the radiative loss function with respect to density and temperature, respectively, evaluated at the equilibrium state, namely

$$\begin{array}{cc}\hfill L_{\rho }& \equiv {\left(\frac{\partial L}{\partial \rho }\right)}_{{\rho }_0,T_0}={\chi }^{\ast }{T}_0^{\alpha }\hfill \end{array}$$(72)

$$\begin{array}{cc}\hfill L_T& \equiv {\left(\frac{\partial L}{\partial T}\right)}_{{\rho }_0,T_0}=\alpha {\rho }_0{\chi }^{\ast }{T}_0^{\alpha -1}.\hfill \end{array}$$(73)

In addition, p′, ρ′, and T′ are related through the perturbed equation of state as

$$\begin{array}{c}\hfill \frac{p\prime }{p_0}=\frac{\rho \prime }{{\rho }_0}+\frac{T\prime }{T_0},\end{array}$$(74)

where the term of order ϵ⁴ has been neglected. In writing this last equation, we assumed, for the sake of simplicity, that the ionization degree of the plasma remains approximately constant. Now, we introduce

$$\begin{array}{c}\hfill {\omega }_{\rho }\equiv \frac{{\rho }_0^2}{p_0}{L}_{\rho },{\omega }_T\equiv \frac{{\rho }_0{T}_0}{p_0}{L}_T.\end{array}$$(75)

In order to obtain the dimensionless energy equation, we introduce the following dimensionless quantities:

$$\begin{array}{cc}& {\overline{k}}_{\Vert ,0}=\frac{k_{\Vert ,0}{T}_0}{c_{\mathrm{A},0}{p}_{0}L},{\overline{c}}_{s,0}=\frac{c_{\mathrm{s},0}}{c_{\mathrm{A},0}},{\overline{\omega }}_{\rho }={\omega }_{\rho }\tau ,{\overline{\omega }}_T={\omega }_T\tau ,\hfill \\ \hfill & \overline{p}\prime =\frac{p\prime }{B_0^2/\mu },\overline{T}\prime =\frac{T\prime }{T_0}\hfill \end{array}$$(76)

Then, the dimensionless second-order equations become

$$\begin{array}{cc}\hfill \frac{\partial \rho \prime }{\partial t}& =-\frac{\partial v_z\prime }{\partial z},\hfill \end{array}$$(77)

$$\begin{array}{cc}\hfill \frac{\partial v_z\prime }{\partial t}& =-\frac{\partial p\prime }{\partial z}-\frac{1}{2}\frac{\partial {B_y\prime }^2}{\partial z},\hfill \end{array}$$(78)

$$\begin{array}{cc}\hfill \frac{\partial p\prime }{\partial t}& =-c_{\mathrm{s},0}^2\frac{\partial v_z\prime }{\partial z}+(\gamma -1)\left[{\kappa }_{\Vert ,0}\frac{{\partial }^2(p\prime -\frac{c_{\mathrm{s},0}^2}{\gamma }\rho \prime )}{\partial z^2}\right]\hfill \\ \hfill & -(\gamma -1)\left[\frac{c_{\mathrm{s},0}^2}{\gamma }({\omega }_{\rho }\rho \prime +{\omega }_T(\frac{\gamma }{c_{\mathrm{s},0}^2}p\prime -\rho \prime ))\right]\hfill \\ \hfill & +(\gamma -1){\eta }_{\mathrm{C},0}{\left(\frac{\partial B_y\prime }{\partial z}\right)}^2,\hfill \end{array}$$(79)

where αω_ρ = ω_T and, using Eqs. (74), (48), and (76), T′ has been written as

$$\begin{array}{c}\hfill T\prime (z,t)=\frac{\gamma }{c_{\mathrm{s},0}^2}p\prime (z,t)-\rho \prime (z,t),\end{array}$$(80)

which has been substituted in Eq. (79). The right-hand side of Eq. (78) involves a quadratic term in $B_y^{\prime }$, which drives slow magnetoacoustic waves, since the second-order magnetic pressure perturbation associated with the Alfvén waves couples with the thermal pressure perturbation associated with the slow waves. In absence of this term and the last term in Eq. (79), Eqs. (77)–(79) would represent free slow magnetoacoustic waves damped by thermal conduction and radiative losses. Here we are not interested in the free slow waves, which have been described elsewhere (Goossens 2003; Roberts 2004, 2006; Priest 2014). Instead, we focus on the solution that represents the generation of slow magnetoacoustic waves due to the nonlinear coupling with the Alfvén waves. On the other hand, from here on it is assumed that $c_{\text{s},\text{0}}{}^2<c_{\text{A},\text{0}}{}^2$; that is to say, we consider a low-β plasma, where the plasma β is defined as

$$\begin{array}{c}\hfill \beta \equiv \frac{c_{\mathrm{s},0}^2}{c_{\mathrm{A},0}^2}=\frac{\gamma p_0}{B_0^2/\mu }·\end{array}$$(81)

In the following subsections, we consider different cases in order to understand how the involved physical mechanisms affect the temporal behavior of second-order perturbations.

3.2.1. Case where radiation and thermal conduction are neglected

First, we consider the case in which radiation and thermal conduction are neglected. We assume an equilibrium temperature of T = 12 000 K and, following Sect. 2.1, the ionization degree in this case is i = 0.82, the mean molecular weight is $\stackrel{\sim }{\mu }=0.73$, and β = 0.07. From Eqs. (77)–(79), the system of equations to be solved is

$$\begin{array}{cc}& \frac{\partial \rho \prime }{\partial t}=-\frac{\partial v_z\prime }{\partial z},\hfill \end{array}$$(82)

$$\begin{array}{cc}& \frac{\partial v_z\prime }{\partial t}=-\frac{\partial p\prime }{\partial z}-\frac{1}{2}\frac{\partial {B_y\prime }^2}{\partial z},\hfill \end{array}$$(83)

$$\begin{array}{cc}& \frac{\partial p\prime }{\partial t}=-c_{\mathrm{s},0}^2\frac{\partial v_z\prime }{\partial z}+(\gamma -1){\eta }_{\mathrm{C},0}{\left(\frac{\partial B_y\prime }{\partial z}\right)}^2.\hfill \end{array}$$(84)

In this system of equations, the only dissipative term is the second term on the right-hand side of Eq. (84), where the expression for $B_y^{\prime }$ obtained from Eq. (70) has been used. The system of equations was solved, numerically, together with the following initial and boundary conditions, respectively,

$$\begin{array}{c}\hfill \rho \prime (z,0)=v_z\prime (z,0)=p\prime (z,0)=0,v_z\prime (z=1,t)=v_z\prime (z=-1,t)=0.\end{array}$$(85)

The case without ambipolar diffusion was studied by Rankin et al. (1994) and Tikhonchuk et al. (1995), among others. To recover it we can neglect the dissipative term in Eq. (84) and combine the remaining equations to obtain a partial differential equation for the longitudinal velocity perturbation, $v_z^{\prime }$,

$$\begin{array}{c}\hfill \frac{{\partial }^2{v}_z\prime }{\partial t^2}-c_{\mathrm{s},0}^2\frac{{\partial }^2{v}_z\prime }{\partial z^2}=-\frac{1}{2}\frac{\partial }{\partial t}\left(\frac{\partial B_y^{\prime 2}}{\partial z}\right).\end{array}$$(86)

Using the D’Alembert’s method, an analytical solution to this equation can be found, namely

$$\begin{array}{c}\hfill v_z\prime =\frac{v_0^2{\omega }_r^{3}sin(2k_{z}z)[{\omega }_{r}sin(2c_{\mathrm{s},0}{k}_{z}t)-c_{\mathrm{s},0}{k}_{z}sin(2{\omega }_{r}t)]}{(8c_{\mathrm{s},0}^3{k}_z^4-8c_{\mathrm{s},0}{k}_z^2{\omega }_r^2)}·\end{array}$$(87)

This analytical solution can be decomposed in two sinusoidal terms with angular frequencies 2c_s,0k_z and 2ω_r, respectively; the longitudinal velocity perturbation comes from the composition of the two sinusoidal curves with different amplitudes and angular frequencies. Considering the low β regime, c_s,0k_z ≪ ω_r, the second term in Eq. (87) can thus be neglected and the dominant angular frequency of the solution is 2c_s,0k_z (Rankin et al. 1994; Tikhonchuk et al. 1995). Using the already known value for c_s,0, the oscillatory period is 132 s, although this value is an approximation to the true period.

To explore the effect of ambipolar diffusion produced by a change of the Alfvén wave wavelength, we consider the following values for L: 10⁷ m, 10⁵ m, and 5 × 10⁴ m; therefore, the corresponding wavelengths are 4 × 10⁷ m, 4 × 10⁵ m, and 2 × 10⁵ m. Using these values for L, three different values for the dimensionless Cowling’s coefficient are obtained, 0.00023, 0.023, and 0.045, respectively, which have been used when numerically solving Eqs. (82)–(84).

Figure 3 shows the temporal and spatial behavior of the heat released due to ambipolar diffusion. We observe that most of the heat is released at the central part of the spatial domain and decreases towards the boundaries, and that it is attenuated in time, as can be expected. When different values of Cowling’s diffusivity are considered, we find that the amount of heat released increases when Cowling’s diffusivity is increased.

Fig. 3. Dimensionless heating versus dimensionless time and space for ηC,0 = 0.023.

Figure 4 displays the temporal behavior of pressure perturbation at z = 0 (left panel) and z = 1 (right panel) for the different values of Cowling’s diffusivity that are considered. At the center (z = 0) of the spatial domain, pressure perturbations are always positive and the amplitude increases when Cowling’s diffusivity is increased. However, for the smallest considered value of Cowling’s diffusivity, pressure perturbations become very slightly negative at the minimum of the oscillations. It seems that the net effect of an increase in Cowling’s diffusivity is that pressure perturbation curves within the positive zone of the plot are slightly shifted up, while a decrease produces a downward shift of the curves, which slightly enter the negative zone of the plot, as can be observed in Fig. 4 (left panel). At the boundary (z = 1) of the spatial domain, we observe two different behaviors: for the highest and the intermediate values of Cowling’s diffusivity, the amplitude of the pressure perturbation goes from negative to positive values, while for the smallest value of Cowling’s diffusivity, the amplitude of the pressure perturbation is negative. This behavior of pressure perturbations is directly related to the behavior of density and temperature perturbations. In the case of density perturbations (Fig. 5), they become negative at the center (z = 0) of the spatial domain for any significant value of Cowling’s diffusivity, while they become positive at the boundary (z = 1). A similar behavior is observed in Fig. 6 for temperature perturbations; we note the sudden increase of the temperature perturbation during the time in which most of the heating is released (Fig. 6, left panel). We also observe that, as expected, as the value of Cowling’s diffusivity increases, so does the temperature perturbation because of the more intense heating. Once the initial temperature increase is over, the temperature perturbations become oscillatory in the same way the pressure and density perturbations do. On the contrary, at z = 1 (Fig. 6, right panel) we only observe oscillatory perturbations without, as before, any sudden increase in temperature due to the lack of heat released at the boundaries of the spatial domain. During the oscillatory behavior, whose amplitude becomes constant after some time, the observed temperature variations are due to the interplay between pressure and density perturbations shown in Eq. (80). Finally, Fig. 7 (left panel) displays the temporal behavior of longitudinal velocity perturbations, at z = 0.5, for the three different values of Cowling’s diffusivity. It is interesting to note that the high-frequency distortions of the sinusoidal shape that appear for the very low value of η_C,0 disappear when the value of η_C,0 is increased. Therefore, a sufficiently intense Cowling’s diffusion is able to damp the high-frequency modulation of the longitudinal velocity perturbation, which is associated with the damped Alfvén wave. Only the low frequency, associated with the slow wave, remains.

Fig. 4. Dimensionless pressure perturbations versus dimensionless time computed at z = 0 (left) and z = 1 (right), for ηC,0 = 0.00023 (black), 0.023 (blue), and 0.045 (red).

Fig. 5. Dimensionless density perturbations versus dimensionless time computed at z = 0 (left) and z = 1 (right), for ηC,0 = 0.00023 (black), 0.023 (blue), and 0.045 (red).

Fig. 6. Dimensionless temperature perturbations versus dimensionless time computed at z = 0 (left) and z = 1 (right), for ηC,0 = 0.00023 (black), 0.023 (blue), and 0.045 (red).

Fig. 7. Left panel: dimensionless longitudinal velocity perturbations versus dimensionless time computed for ηC,0 = 0.00023 (black), 0.023 (blue), and 0.045 (red); Right panel: dimensionless longitudinal velocity perturbation versus dimensionless time. Full analytical solution obtained from Eq. (89) (red) and numerical solution (blue) (z = 0.5, v0 = 0.15, kz = π/2, ηC,0 = 0.045 and 0.00023).

The temporal behavior of pressure, density, and temperature perturbations can be understood as follows. For any significant value of Cowling’s diffusivity, the release of heat due to ambipolar diffusion produces, for short periods, an increment of temperature around the center of the spatial domain. Later on, the behavior of the temperature becomes oscillatory. The rise in temperature at the center of the spatial domain, because of the heating, causes a local increment of pressure. In turn, this over-pressure pushes the material away from the center, so that a longitudinal flow towards the boundaries is generated. As a result of this, mass density decreases at the center and increases near the boundaries (see Fig. 8). We note that this effect is the opposite of that of the ponderomotive force, which tends to accumulate mass at the center of the domain. Which of the two effects dominates depends on the efficiency of the heating due to ambipolar diffusion.

Fig. 8. Dimensionless density perturbation versus dimensionless time and space for ηC,0 = 0.045.

On the other hand, Eqs. (82)–(84) can be combined to obtain the following wave equations:

$$\begin{array}{cc}\hfill \frac{{\partial }^{2}p\prime }{\partial t^2}-c_{\mathrm{s},0}^2\frac{{\partial }^{2}p\prime }{\partial z^2}& =\frac{c_{\mathrm{s},0}^2}{2}\frac{\partial }{\partial z}\left(\frac{\partial B_y^{\prime 2}}{\partial z}\right)\hfill \\ \hfill & +(\gamma -1){\eta }_{\mathrm{C},0}\frac{\partial }{\partial t}{\left(\frac{\partial B_y\prime }{\partial z}\right)}^2\hfill \end{array}$$(88)

$$\begin{array}{cc}& \frac{{\partial }^2{v}_z\prime }{\partial t^2}-c_{\mathrm{s},0}^2\frac{{\partial }^2{v}_z\prime }{\partial z^2}=-\frac{1}{2}\frac{\partial }{\partial t}\left(\frac{\partial B_y^{\prime 2}}{\partial z}\right)\hfill \\ \hfill & -(\gamma -1){\eta }_{\mathrm{C},0}\frac{\partial }{\partial z}{\left(\frac{\partial B_y\prime }{\partial z}\right)}^2.\hfill \end{array}$$(89)

Then, analytical solutions for pressure (p′) and longitudinal velocity perturbations ($v_z^{\prime }$) can be found by means of D’Alembert’s method, and for density perturbations (ρ′) from Eq. (82). These solutions are much more complex than in the case without ambipolar diffusion (Eq. (87)) and are not shown here. Instead, in Fig. 7 (right panel), we compare the obtained analytical solution for velocity perturbation with the previously obtained numerical solution for two different values of Cowling’s diffusivity and find a perfect agreement.

These full analytical solutions for the pressure, longitudinal velocity, and density perturbations obtained from Eqs. (88), (89), and (82) can be used to compute the temporal behavior of the internal energy, as well as the energy associated with Alfvén and slow waves. In dimensionless form, the change in the plasma’s total internal energy is computed from

$$\begin{array}{c}\hfill \mathrm{\Delta }{E}_{\mathrm{int}}={\displaystyle {\int }_{-1}^1\frac{p(z)-p_0}{\gamma -1}\mathrm{d}z,}\end{array}$$(90)

where p(z)−p₀ = ϵ²p′ is the final pressure distribution and p₀ is the initial constant pressure. The change in the slow wave energy is

$$\begin{array}{c}\hfill \mathrm{\Delta }{E}_{\mathrm{slow}}=\frac{1}{2}{\displaystyle {\int }_{-1}^1(v_z^{\prime 2}+\frac{p^{\prime 2}}{c_{\mathrm{s},0}^2})\mathrm{d}z,}\end{array}$$(91)

which involves kinetic and thermal energy; for Alfvén waves it is

$$\begin{array}{c}\hfill \mathrm{\Delta }{E}_{\mathrm{Alfven}}=\frac{1}{2}{\displaystyle {\int }_{-1}^1(v_y^{\prime 2}+B_y^{\prime 2})\mathrm{d}z,}\end{array}$$(92)

which involves kinetic and magnetic energy. Figure 9 shows the temporal behavior of the variation of these energies. It is clear that almost all the energy from Alfvén waves goes to increment the internal energy and, as a consequence, to increment the plasma temperature. On the contrary, only a negligible amount of energy goes to slow waves. In order to obtain an approximation to the plasma’s final temperature, we neglect the amount of energy transferred from Alfvén to slow waves and we assume

$$\begin{array}{c}\hfill E_{\mathrm{KA}}=E_{\mathrm{int}},\end{array}$$(93)

Fig. 9. Dimensionless variation of the waves and internal energies versus dimensionless time. Black: internal energy; red: Alfvén wave energy; blue: slow wave energy. (T = 12 000 K, v0 = 0.15, kz = π/2, ηC,0 = 0.0023).

where E_KA is the total kinetic energy of the Alfvén wave given by

$$\begin{array}{c}\hfill E_{\mathrm{KA}}={\displaystyle {\int }_{-1}^1\frac{1}{2}{\rho }_0{v}_y^{\prime 2}\mathrm{d}z=\frac{1}{2}{\rho }_0{v}_0^2,}\end{array}$$(94)

where ${\mathit{v}}_{\mathit{y}}\prime ={\mathit{v}}_{0}cos(\frac{\pi }{2}z),$ the initial perturbation for the Alfvén wave. Therefore, the above expression provides us with the initial energy in the Alfvén wave. The E_int is the internal energy given by

$$\begin{array}{c}\hfill E_{\mathrm{int}}={\displaystyle {\int }_{-1}^1\frac{p(z)-p_0}{\gamma -1}\mathrm{d}z=2\frac{p-p_0}{\gamma -1},}\end{array}$$(95)

where p(z) is the final pressure distribution and which, for the sake of simplicity, has been assumed as uniform and equal to p. Then, substituting the energy expressions in Eq. (93), we can obtain an expression for the final temperature, namely:

$$\begin{array}{c}\hfill T=T_0(1+\frac{1}{4}\gamma (\gamma -1)\frac{v_0^2}{c_{\mathrm{s},0}^2}).\end{array}$$(96)

Then,

$$\begin{array}{c}\hfill \mathrm{\Delta }T=T-T_0=T_0\left(\frac{1}{4}\gamma (\gamma -1)\frac{v_0^2}{c_{\mathrm{s},0}^2}\right),\end{array}$$(97)

dividing by T₀ and introducing c_A,0, the temperature variation can be written in dimensionless form as

$$\begin{array}{c}\hfill \mathrm{\Delta }T=\frac{1}{4}\gamma (\gamma -1)\frac{v_0^2}{c_{\mathrm{s},0}^2}·\end{array}$$(98)

Using T₀ = 12 000 K for the equilibrium temperature, η_C,0 = 0.045, and the equilibrium Alfvén and sound speeds, we obtain ΔT = 0.108. On the other hand, using the above-mentioned analytical solutions for pressure and density perturbations, Eq. (80), and the same physical conditions, we can compute the spatial average of temperature perturbation by means of

$$\begin{array}{c}\hfill ⟨T\prime [z,t]⟩=\frac{1}{2L}{\displaystyle {\int }_{-L}^{L}T\prime [z,t]\mathrm{d}z.}\end{array}$$(99)

This computation is done for a time at which all the Alfvén wave’s energy has been converted to internal energy and the obtained result fully coincides with the value of ΔT that was obtained before.

3.2.2. Effect of radiation

In this case, our energy equation comes from Eq. (79), when thermal conduction is neglected and a radiative loss function (Hildner 1974) is considered. Then, the system of equations to be solved is:

$$\begin{array}{cc}\hfill \frac{\partial \rho \prime }{\partial t}& =-\frac{\partial v_z\prime }{\partial z},\hfill \end{array}$$(100)

$$\begin{array}{cc}\hfill \frac{\partial v_z\prime }{\partial t}& =-\frac{\partial p\prime }{\partial z}-\frac{1}{2}\frac{\partial {B_y\prime }^2}{\partial z},\hfill \end{array}$$(101)

$$\begin{array}{cc}\hfill \frac{\partial p\prime }{\partial t}& =-c_{\mathrm{s},0}^2\frac{\partial v_z\prime }{\partial z}+(\gamma -1)[-\frac{c_{\mathrm{s},0}^2}{\gamma }({\omega }_{\rho }\rho \prime +{\omega }_T(\frac{\gamma }{c_{\mathrm{s},0}^2}p\prime -\rho \prime ))]\hfill \\ \hfill & +(\gamma -1)\left[{\eta }_{\mathrm{C},0}{\left(\frac{\partial B_y\prime }{\partial z}\right)}^2\right]·\hfill \end{array}$$(102)

This system was solved together with the same initial and boundary conditions as in Sect. 3.2.1. In order to determine the effect produced by radiative losses, we consider three different equilibrium temperatures, T = 12 000 K, T = 8000 K, and T = 6000 K, typical of quiescent solar prominences, together with two different values for L, 10⁵ and 10⁷ m, in an attempt to identify the effect of a change in the wavelength, λ = 4L, of the Alfvén wave. Next, the dimensionless numerical values of Cowling’s diffusity, η_C,0, are computed together with ω_T and ω_ρ. Figures 10–11 show the temporal behavior of density, pressure, longitudinal velocity, and temperature perturbations for L = 10⁵ m. The results show that for L = 10⁵ m, the shorter wavelength, the perturbations corresponding to T = 12 000 K are quickly damped; for L = 10⁷ m, the longer wavelength, the opposite happens and the perturbations corresponding to T = 6000 K are more quickly damped. This behavior can be understood following Soler (2010). In this work, the author obtained an analytical approximation for the imaginary part of the angular frequency, in the case of nonadiabatic slow waves in a filament thread. When radiation and thermal conduction are considered, the imaginary part of the angular frequency is given by

$$\begin{array}{c}\hfill {\omega }_I=-\frac{\gamma -1}{2\gamma }\frac{(\gamma -1)({\kappa }_{\Vert ,0}{k}_z^2+{\omega }_T)+{\omega }_{\rho }}{c_{\mathrm{s},0}^2{k}_z^2+{(\gamma -1)}^2{({\kappa }_{\Vert ,0}{k}_z^2+{\omega }_T)}^2}{c}_{\mathrm{s},0}^2{k}_z^2.\end{array}$$(103)

Fig. 10. Left panel: dimensionless density perturbations versus dimensionless time, with ambipolar diffusion and radiative losses included. Right panel: dimensionless pressure perturbations versus dimensionless time, with ambipolar diffusion and radiative losses included (z = 0, L = 105 m, T = 6000 K and ηC,0 = 0.03 (black curve), 8000 K and ηC,0 = 0.026 (blue curve), 12 000 K and ηC,0 = 0.023 (red curve)).

Fig. 11. Left panel: dimensionless longitudinal velocity perturbations versus dimensionless time, with ambipolar diffusion and radiative losses included computed at z = 0.5. Right panel: dimensionless temperature perturbations versus dimensionless time, with ambipolar diffusion and radiative losses included computed at z = 0 (L = 105 m, T = 6000 K and ηC,0 = 0.03 (black curve), 8000 K and ηC,0 = 0.026 (blue curve), 12 000 K and ηC,0 = 0.023 (red curve)).

Then, once this quantity has been computed for the different temperatures, the damping time can be computed from

$$\begin{array}{c}\hfill {\tau }_D=\frac{1}{|{\omega }_I|}·\end{array}$$(104)

When only radiation is considered (κ_∥ = 0), with the help of the Eqs. (103) and (104), we can obtain the damping times for each of the considered equilibrium temperatures and different values of L (see Table 2). For L = 10⁷ m, the results show that the perturbations corresponding to the lower temperature are the ones that are more quickly damped. Conversely, when L is decreased to 10⁵ m, the perturbations corresponding to the higher temperature are the ones that are more quickly damped. When only radiation is considered and the temperature is kept fixed, all the quantities in Eq. (103) are constants except the wavenumber, k_z, which depends on L. Then, Eq. (104) can be written as

$$\begin{array}{c}\hfill {\tau }_D=F+G{L}^2,\end{array}$$(105)

suggesting a parabolic behavior of the damping time versus L, while F and G are constants. This can be seen in Fig. 12: for large values of L, no crosses between the different curves appear and the damping time is always much longer for high temperatures than for low ones. However, for small values of L, we observe the presence of crosses between the different curves and, as a consequence, the damping time for high temperatures is shorter than for low temperatures. For a constant temperature, the behavior of the damping time only depends on L (see Table 2), a quantity which is related to the longitudinal wavenumber through $k_z=\frac{\pi }{2L}$.

Fig. 12. Logarithm of the damping time versus L for T = 12 000 K (red), T = 8000 K (blue), T = 6000 K (black).

On the other hand, for L = 10⁵ m, the timescale is τ = 1.58 s, which means that the dimensionless time interval considered in Figs. 10 and 11 corresponds to 94.8 s. Then we can observe that: The perturbations are completely damped for a temperature of 12 000 K with a damping time of 49 s; they are slightly damped for 8000 K with a damping time of 104 s; and they are nearly undamped for 6000 K, with a damping time of 627 s. When L = 10⁷ m is considered, the behavior of the temporal perturbation damping is reversed. The comparison between the obtained damping times (Table 2) and the temporal behavior of perturbations in Figs. 10 and 11 confirms that the analytical approximation is valid and explains the temporal attenuation of perturbations when different equilibrium temperatures, as well as different characteristic length scales, L, are considered.

Furthermore, it is also interesting to compare the damping times for Alfvén and slow waves. For L = 10⁷ m, the damping times of Alfvén waves are much greater than those of slow waves, while the opposite happens when L = 10⁵ m. Modifying the value of L means modifying the longitudinal wavenumber k_z; therefore, when L is increased the wavenumber decreases and, following Eq. (68), the damping time is increased. Conversely, when L is decreased the damping time decreases. These results show that when the wavelength of the Alfvén wave, 4L, is increased, the efficiency of the ambipolar diffusion decreases and the opposite happens when is decreased. Then, when ambipolar diffusion is very efficient, Alfvén waves are quickly damped while slow waves still remain in the plasma; however, when ambipolar diffusion is inefficient, slow waves are more quickly damped. Figure 13 displays the temporal behavior of Alfvén, internal, and slow wave energies, showing how these energies are damped in time due to ambipolar and radiative losses. Then, taking into account the temporal behavior of the energies, we can explain the temporal behavior of the perturbations. At the beginning, they should behave as in the case without radiation; however, once radiative losses become dominant, these perturbations must be damped in time. Finally, Fig. 14 (left panel) shows the spatial distribution of the two terms contributing to the energy balance separately and for a fixed time; it can be seen that, for the chosen time, the injected heating is still dominant and radiative losses are small. Therefore, the balance is positive and the plasma is heated.

Fig. 13. Dimensionless variation of internal (red), slow (blue) and Alfvén (black) wave energies versus dimensionless time (T = 6000 K and ηC,0 = 0.003).

3.2.3. Combined effect of radiation and thermal conduction

In this case, we must solve Eqs. (77)–(79) together with the following initial and boundary conditions:

$$\begin{array}{cc}& \rho \prime (z,0)=v_z\prime (z,0)=p\prime (z,0)=0,\hfill \\ \hfill & v_z\prime (z=1,t)=v_z\prime (z=-1,t)=0,\hfill \\ \hfill & {\partial }_{z}p\prime (z=1,t)={\partial }_{z}p\prime (z=-1,t)=0,\hfill \\ \hfill & {\partial }_z\rho \prime (z=1,t)={\partial }_z\rho \prime (z=-1,t)=0.\hfill \end{array}$$(106)

Figure 14 (right panel) shows, at a fixed time, the spatial distribution of the balance between contributions coming from ambipolar diffusion, radiation, and thermal conduction; we observe that radiation dominates the energy loss while thermal conduction losses are negligible. This conclusion is in agreement with the temporal behavior of the different perturbations, which is quite similar to the case when only radiative losses are considered.

Fig. 14. Left panel: spatial distribution of the energy balance. Red curve: radiative losses; black curve: heating due to ambipolar diffusion; blue curve: radiative losses plus heating due to ambipolar diffusion (t = 1, T = 6000 K and ηC,0 = 0.003). Right panel: spatial distribution of the energy balance. Red curve: radiative losses; blue curve: thermal conduction losses; black curve: heating due to ambipolar diffusion; brown curve: radiative plus thermal conduction losses plus heating due to ambipolar diffusion (t = 3, T = 6000 K and ηC,0 = 0.0023).

3.3. Third-order equations

The third-order equations represent nonlinear corrections to the Alfvén modes. The dimensionless third-order equations are:

$$\begin{array}{cc}& \frac{\partial v_{3y}\prime }{\partial t}=-v_z\prime \frac{\partial v_y\prime }{\partial z}+\frac{\partial B_{3y}\prime }{\partial z}-\rho \prime \frac{\partial v_y\prime }{\partial z},\hfill \end{array}$$(107)

$$\begin{array}{cc}& \frac{\partial B_{3y}\prime }{\partial t}=\frac{\partial v_{3y}\prime }{\partial z}-\frac{\partial (B_y\prime v_z\prime )}{\partial z}+{\eta }_{\mathrm{C},0}\frac{{\partial }^2{B}_{3y}\prime }{\partial z^2}·\hfill \end{array}$$(108)

These can be combined to obtain a wave equation for $v_{3y}^{\prime }$ so that