© The Authors 2024

1. Introduction

Solar quiescent prominences are majestic structures of plasma sustained against gravity thanks to magnetic fields, whose strengths vary between 3 G and 60 G (see, e.g., Mackay et al. 2010; Gibson 2018). Such structures can be seen as bright structures in the solar limb or as dark structures on the solar disk, in the case of which they are referred to as filaments that are formed in filament channels. Projected on the photosphere, solar prominences are near to a neutral line, which separates regions of opposite magnetic polarity (Babcock & Babcock 1955; Howard & Harvey 1964). Typically, solar prominences have temperatures between 7500 K and 9000 K, gas pressure between 0.1 dyn cm⁻² and 1 dyn cm⁻², and number densities between 10¹⁰ cm⁻³ and 10¹¹ cm⁻³ (Labrosse et al. 2010; Parenti 2014; Engvold 2015). Moreover, solar prominences are usually suspended at heights that corresponds to the lower corona. In such environment, temperatures are above 10⁶ K, but the number density is low, ∼10⁸ cm⁻³ (see, e.g., Priest 2014). Consequently, the plasma inside a prominence is cold and dense compared with the coronal plasma. Importantly, the plasma in prominences is only partially ionized (Labrosse et al. 2010).

High-cadence and high-resolution Hα observations showed that prominences are formed by a myriad of substructures called prominence threads (Lin et al. 2005, 2007, 2008). Prominence threads are quite thin (∼0.21 Mm) and long (between ∼3.5 and 28 Mm) structures, which are believed to be embedded in much longer magnetic tubes (see, e.g., Lin et al. 2005, 2007) of lengths of the order of 100 Mm or more whose feet are rooted in the lower atmosphere (Terradas et al. 2008a). Consequently, prominence threads are usually modeled as thin magnetic flux tubes with their ends fixed at the photosphere (Ballester & Priest 1989; Rempel et al. 1999; Soler et al. 2012; Adrover-González et al. 2021; Melis et al. 2023).

Prominence threads are so dynamic that their continuous presence in Hα observations is typically short. It ranges between few minutes and 20 min (Lin et al. 2005). Magnetohydrodynamic (MHD) waves and oscillations play a predominant role in the dynamics of threads (see Arregui et al. 2018). There are many reports of transverse oscillations and propagating waves in prominence threads with periods between 2 min and ∼15 min (Lin et al. 2007, 2009; Ning et al. 2009; Orozco Suárez et al. 2014; Li et al. 2022). These oscillations are interpreted as kink MHD waves (Terradas et al. 2008a; Lin et al. 2009; Soler et al. 2011; Okamoto et al. 2015; Li et al. 2022). The driver of these waves is believed to be the convective motions at the photosphere, where the prominence magnetic field is anchored (Hillier et al. 2013).

The present paper deals with torsional Alfvén waves, which are transverse MHD waves that periodically twist the magnetic field lines in a flux tube. It has been postulated that they might play a role in the heating of the solar corona (Hollweg 1978; Cranmer & van Ballegooijen 2005; Cargill & de Moortel 2011; Mathioudakis et al. 2013; Soler et al. 2019; Van Doorsselaere et al. 2020; Nakariakov & Kolotkov 2020) and in the acceleration of the solar wind (Charbonneau & MacGregor 1995; Cranmer 2009; Matsumoto & Suzuki 2012; Shoda et al. 2018). In straight flux tubes with a purely axial magnetic field, these waves can maintain their pure magnetic nature despite inhomogeneities (Goossens et al. 2011). Torsional Alfvén waves are nearly incompressible and are polarized perpendicularly to the magnetic field lines producing periodic perturbations in the azimuthal components of the magnetic field and the velocity. Recently, Soler et al. (2019, 2021) found that torsional Alfvén waves can transport large energy fluxes when they propagate from the photosphere to a coronal loop despite the filtering role of the chromosphere. They found that the energy flux is channeled at the frequencies that match the natural frequencies of the coronal loop, generating global standing torsional oscillations. The generation of standing modes is possible because coronal loops can act as Alfvén cavity resonators (Hollweg 1984a,b). However, the energy distribution is not equally distributed among the different eigenmodes, with the fundamental mode being the major contributor.

To the best of our knowledge, there is no direct evidence of torsional Alfvén waves in prominence threads, although there have been direct reports of rotational motions caused by magnetic reconnection (Okamoto et al. 2016). Nonetheless, Jess et al. (2009) detected torsional Alfvén waves at photospheric bright points through spectral line nonthermal widening (see Zaqarashvili 2003). Morton et al. (2013) reported the presence of torsional Alfvén and kink waves in chromospheric swirls. De Pontieu et al. (2012, 2014a) detected torsional motions in spicules that could be compatible with torsional Alfvén waves. This type of waves were also reported in the chromosphere by Srivastava et al. (2017). More recently, Kohutova et al. (2020) detected torsional Alfvén waves at coronal heights using the Interface Region Imaging Spectrograph (IRIS; De Pontieu et al. 2014b). Moreover, some oscillations generated during solar flares can be interpreted as torsional Alfvén waves (Aschwanden & Wang 2020) and torsional Alfvén waves have been reported in the solar wind (Raghav & Kule 2018). Therefore, in view of their ubiquity, it is likely that torsional Alfvén waves are also present in prominences, although their direct detection remains elusive.

In a line-tied flux tube that is transversely nonuniform, either in density or in magnetic field, there is an Alfvén frequency continuum (see Halberstadt & Goedbloed 1993). Consequently, Alfvén modes at different radial positions in the flux tube become out of phase as time increases. This phenomenon is the well-known process of phase-mixing (see, e.g., Heyvaerts & Priest 1983; Nocera et al. 1984; De Moortel et al. 2000; Smith et al. 2007; Prokopyszyn et al. 2019; Díaz-Suárez & Soler 2021a). Phase-mixing is a linear process that generates shear flows perpendicularly to the magnetic field lines and transports Alfvén wave energy from large to small perpendicular scales. The rhythm at which the energy is transported depends on the gradient of the Alfvén frequency (see Mann et al. 1995). Eventually, the phase-mixing shear flows trigger the Kelvin–Helmholtz instability (KHI) as Heyvaerts & Priest (1983) and Browning & Priest (1984) analytically predicted. This result is numerically confirmed in Guo et al. (2019) and Díaz-Suárez & Soler (2021b). Although with some important differences, an equivalent dynamics also appears in simulations of kink waves (see, e.g., Terradas et al. 2008b, 2018; Antolin et al. 2014; Magyar & Van Doorsselaere 2016; Howson et al. 2017a,b; Karampelas et al. 2017, 2019; Karampelas & Van Doorsselaere 2018; Guo et al. 2019; Pascoe et al. 2020; Shi et al. 2021; Magyar et al. 2022). The KHI generates vortices that nonlinearly break into smaller and smaller vortices leading naturally to turbulence. Although not related to the Alfvén waves, there are observations of the KHI in coronal mass ejections (Foullon et al. 2011) and in prominences (Berger et al. 2017; Hillier & Polito 2018; Yang et al. 2018). In the latter case, turbulence is also reported (see Leonardis et al. 2012).

In the present work we study the nonlinear evolution of torsional Alfvén waves in quiescent prominence threads. The purpose of this investigation is twofold. On the one hand, we aim to the explore the process of turbulence generation mediated by torsional waves. Such a mechanism has been studied before in the case of coronal loops (Díaz-Suárez & Soler 2021b, 2022) but, to the best of our knowledge, not in prominence threads. Unlike the fully ionized coronal loops, prominence threads are only partially ionized. In the partially ionized prominence plasma, ambipolar diffusion and, to a lesser extent, Ohmic diffusion, are important dissipation mechanisms (Khomenko et al. 2014a; Ballester et al. 2018; Melis et al. 2021). Consequently, ambipolar and Ohmic diffusion are included here as nonideal effects that can dissipate Alfvén waves and the associated turbulence, and so potentially heat prominence threads. On the other hand, we aim to explore what kind of signatures the dynamics of torsional oscillations may leave in synthetic Hα observations, so that observers may look for those signatures in real observations in order to detect the still elusive torsional waves in prominences.

This paper is organized as follows. In Sect. 2, we describe the numerical setup. The results from both ideal MHD and nonideal MHD simulations are given in Sect. 3. The synthetic Hα observations are presented and analyzed in Sect. 4. Finally, we discuss the conclusions of this work in Sect. 5.

2. Numerical setup

2.1. Prominence thread model

We represent a prominence thread as a straight magnetic flux tube of length, L, and radius, R, permeated by a uniform longitudinal magnetic field, namely, $\mathit{B}=B_0\widehat{z}$, where B₀ = 5.2 G throughout. We consider R = 1 Mm and L = 50R, so that L = 50 Mm. The value of L used here is shorter than that usually reported from observations by a factor of 2, approximately. Terradas et al. (2008a) inferred the minimum length of the magnetic tube of an active-region prominence thread reported by Okamoto et al. (2007) as L ≈ 104.8 Mm. The reason for considering a shorter tube is to speed up the computations, since the periods of the standing waves are proportional to L.

The cool and dense plasma of the prominence thread is surrounded by the hot and light coronal plasma. For simplicity, we ignore gravity and specify the equilibrium density, ρ₀, that varies in both the radial, r, and longitudinal, z, directions namely,

$$\begin{array}{c}\hfill {\rho }_0(r,z)=\{\begin{array}{ccc}{\rho }_{\mathrm{i}}(z),\hfill & \mathrm{if}\hfill & r\le R-\frac{l}{2},\hfill \\ {\rho }_{\mathrm{tr}}(r,z),\hfill & \mathrm{if}\hfill & R-\frac{l}{2}<r<R+\frac{l}{2},\hfill \\ {\rho }_{\mathrm{e}},\hfill & \mathrm{if}\hfill & r\ge R+\frac{l}{2},\hfill \end{array}\end{array}$$(1)

where ρ_i(z) is the internal density that varies along the tube, ρ_e is the external coronal density assumed uniform, and ρ_tr(r, z) is the density in a transversely nonuniform layer of thickness, l, that connects the internal and external plasmas. In this work we used l/R = 0.6. After Soler et al. (2015), we adopt a Lorentzian profile along the flux tube for the internal density, namely,

$$\begin{array}{c}\hfill {\rho }_{\mathrm{i}}(z)=\frac{{\rho }_{\mathrm{i},0}}{1+4(\chi -1)z^2/L^2}·\end{array}$$(2)

In Eq. (2), ρ_i, 0 is the density at the centre of the flux tube, z = 0, and χ ≡ ρ_i, 0/ρ_i(z = L/2) is the ratio of the central density to the footpoint density. The larger the value of χ, the more concentrated the distribution of the density is around the center of the flux tube (see, e.g., Fig. 2 of Martínez-Gómez et al. 2022). In turn, the density in the transversely nonuniform layer is prescribed as

$$\begin{array}{c}\hfill {\rho }_{\mathrm{tr}}(r,z)=\frac{{\rho }_{\mathrm{i}}(z)}{2}\{[1+\frac{{\rho }_{\mathrm{e}}}{{\rho }_{\mathrm{i}}(z)}]-[1-\frac{{\rho }_{\mathrm{e}}}{{\rho }_{\mathrm{i}}(z)}]sin\left[\frac{\pi }{l}(r-R)\right]\}·\end{array}$$(3)

In our background model we used ρ_e = 5.02 × 10⁻¹³ kg m⁻³, ρ_i, 0 = 100ρ_e = 5.02 × 10⁻¹¹ kg m⁻³, and χ = 100. Figure 1 shows a sketch of the prominence thread model and Fig. 2 shows one-dimensional longitudinal and radial cuts of the equilibrium density profile.

Fig. 1. Sketch of the prominence thread model. The solar photosphere is represented by the two gray planes at both ends of the magnetic flux tube.

Fig. 2. Equilibrium density profile. Top panel: longitudinal dependence at r = 0. Bottom panel: transverse dependence at y = z = 0. We used l/R = 0.6, ρi, 0 = 100ρe, L/R = 50, and χ = 100. The density profiles are normalized with respect to the external, coronal density.

The background gas pressure in the model, p₀, is uniform. In the realistic corona and prominences, the plasma is magnetically dominated, so that the plasma $\beta =2p_0{\mu }_0/B_0^2\ll 1$, where μ₀ is the vacuum magnetic permeability. In our model, we chose the value of p₀ so that β = 0.048. The equilibrium Alfvén speed, $v_{\mathrm{A},0}(r,z)=B_0/\sqrt{{\mu }_0{\rho }_0(r,z)}$, and the equilibrium sound speed, $c_{\mathrm{s},0}(r,z)=\sqrt{\gamma p_0/{\rho }_0(r,z)}$, where γ is the adiabatic constant, are displayed in Fig. 3 in longitudinal and radial cuts to the thread. We can see a sharp variation in v_A, 0 and c_s, 0 across the prominence thread owing to the large density contrast between the corona and the core of the thread. The variation of both speeds is smoother in the longitudinal direction. The external value of the Alfvén speed is our reference velocity, namely v_A, e = 654.7 km s⁻¹.

Fig. 3. Same as Fig. 2, but for the sound and Alfvén speeds in the equilibrium. In both panels, the blue solid line corresponds to the Alfvén speed and the red dashed line corresponds to the sound speed. Both speeds are normalized with respect to the external Alfvén speed, vA, e.

2.2. Numerical code

We performed time-dependent numerical simulations with the PLUTO code (Mignone et al. 2007), which solves the 3D resistive MHD equations with a finite-volume, shock-capturing spatial discretization. The equations solved by PLUTO are as follows:

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}=-\mathrm{\nabla }·\left(\rho \mathit{v}\right),\hfill \end{array}$$(4)

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

$$\begin{array}{cc}& \frac{\partial \mathit{B}}{\partial t}=\mathrm{\nabla }\times (\mathit{v}\times \mathit{B})-{\mu }_0\mathrm{\nabla }\times (\widehat{\eta }·\mathit{j}),\hfill \end{array}$$(6)

$$\begin{array}{cc}& \frac{\mathrm{D}p}{\mathrm{D}t}=\frac{\gamma p}{\rho }\frac{\mathrm{D}\rho }{\mathrm{D}t}+(\gamma -1){\mu }_0(\widehat{\eta }·\mathit{j})·\mathit{j}.\hfill \end{array}$$(7)

In Eqs. (4)–(7), $\frac{\mathrm{D}}{\mathrm{D}t}\equiv \frac{\partial }{\partial t}+\mathit{v}·\mathrm{\nabla }$ denotes the total derivative, ρ is the mass density, v is the velocity, B is the magnetic field, p is the gas pressure, and j = (∇×B)/μ₀ is the current density. Moreover, $\widehat{\eta }$ is the resistivity tensor, which is defined in PLUTO as a diagonal tensor in the Cartesian coordinate frame, namely

$$\begin{array}{c}\hfill \widehat{\eta }=\left(\begin{array}{ccc}{\eta }_x& 0& 0\\ 0& {\eta }_y& 0\\ 0& 0& {\eta }_z\end{array}\right),\end{array}$$(8)

with η_x, η_y, and η_z the x-, y-, and z-components of resistivity, respectively.

The PLUTO code solves Eqs. (4)–(7) in Cartesian coordinates. We used a uniform grid of 800 × 800 × 250 points. We performed the simulations in a computational box where x/R ∈ [−3,3], y/R ∈ [−3,3], and z/R ∈ [−25,25]. Therefore, the spatial resolution in x- and y-directions is 7.5 km while in z-direction is 200 km. We used a fifth-order weighted essentially non-oscillatory (WENO) algorithm for spatial reconstruction (Borges et al. 2008) and a Roe-Riemann (Roe 1981) solver to compute the numerical fluxes. Moreover, we used the hyperbolic divergence cleaning technique (Dedner et al. 2002) to maintain the solenoidal constraint on the magnetic field, which couples the divergence of the magnetic field and Eq. (6) to generalized Lagrange multipliers. Because the magnetic field is current-free and force-free, we used the background-field-splitting technique (Powell 1994), which only evolves the magnetic field perturbations. We perform both ideal and resistive MHD simulations. In the ideal simulations ($\widehat{\eta }=0$), an explicit total variation diminishing third-order Runge-Kutta algorithm is used for the temporal evolution. In the resistive simulations ($\widehat{\eta }\ne 0$), the explicit time step becomes too small, as it depends quadratically on the grid size and inversely on the diffusion coefficient. In this case, we use the super-time stepping technique (STS; Alexiades et al. 1996) implemented in PLUTO. STS allows us to save computational time by accelerating the explicit temporal scheme and requiring stabilization only at the advective time super-steps, but not in the various substeps in which each super-step is divided into.

Although the temperature, T, is not a variable directly evolved by the PLUTO code, we are interested in studying its evolution over the span of the simulations. To this end, we implemented in PLUTO the computation of the temperature through an iterative method using the local values of density and gas pressure, so that T is a secondary, user-defined variable of the code. In addition, the prominence plasma is partially ionized (see, e.g., Parenti 2014). Therefore, we also need to determine the plasma ionization fraction, ξ_i, defined as the ratio of the ion density to the total density. At every time step and at each grid point, we solved the equation of state for a partially ionized hydrogen plasma, namely,

$$\begin{array}{c}\hfill p=(1+{\xi }_i)\rho \stackrel{\sim }{R}T,\end{array}$$(9)

where $\stackrel{\sim }{R}$ is the ideal gas constant. In the prominence plasma, ξ_i is a function of the local pressure and temperature. We used the tabulated values of ξ_i given in Table 1 of Heinzel et al. (2015) for a height of 10 Mm over the photosphere. Equation (9) is coupled with the values of the table. Starting with the assumption that ξ_i = 1 (full ionization), Eq. (9), together with the data in Table 1 of Heinzel et al. (2015), is solved iteratively until the computed values of T and ξ_i converge. Since the table only contains limited ranges of pressure and temperature, the following conditions are applied in the process when the values are outside those of the table. Using the results from Gouttebroze & Labrosse (2009), for temperatures higher than 2 × 10⁴ K, we assumed full ionization. If the pressure is larger or smaller than the maximum or minimum pressure value from the table, or if the temperature is smaller than 6000 K, we saturated to the closest value in the table. However, we note that the condition for the pressure was never actually applied in the simulations included here, as the minimum and maximum pressure values found during the simulations fell within the range of the table in all cases. Particularly, the pressure was found to range between 0.046 and 0.057 dyn cm⁻², approximately. Typically, prominences are composed of ∼10% helium. However, for simplicity, we assumed a plasma composed solely of hydrogen.

We used an initial perturbation with the aim of exciting the longitudinally fundamental mode of standing torsional Alfvén waves. As in Díaz-Suárez & Soler (2021b), we perturbed the azimuthal component of velocity, as follows:

$$\begin{array}{c}\hfill \mathit{v}(t=0)=v_{0}A(r)cos\left(\frac{\pi }{L}z\right)\widehat{\phi },\end{array}$$(10)

where v₀ is the maximum velocity amplitude and A(r) contains the radial dependence (see details in Díaz-Suárez & Soler 2021b). A radial cut of the normalized azimuthal component of velocity at the centre of the tube can be seen in top panel of Fig. 4 of Díaz-Suárez & Soler (2022; see the red dot-dashed line). We set v₀ = 0.01 v_A, e, so that v₀ = 6.547 km s⁻¹, which is of the order of the velocity amplitude in small-amplitude prominence oscillations (Lin et al. 2009; Arregui et al. 2018). The longitudinal dependence as $cos\left(\frac{\pi }{L}z\right)$, set in the initial condition of Eq. (10), excites the fundamental mode as well as other longitudinal harmonics with even symmetry with respect to z = 0. These additional harmonics will be present because, in a longitudinally nonuniform tube, the spatial dependence of the eigenmodes deviates from the canonical harmonic dependence (see, e.g., Andries et al. 2005; Soler et al. 2015). However, in the simulations the evolution is largely dominated by the dynamics of the longitudinally fundamental mode.

Regarding the boundary conditions, we used outflow conditions, that is, zero gradient for all the variables in all lateral boundaries. On the top and bottom boundaries, z = ±L/2, we imposed line-tying conditions so as to mimic the anchoring of the field lines in the solar photosphere. To perform this task, we fixed the three components of velocity to zero while the z-component of the magnetic field is fixed to the equilibrium value. The remaining variables, namely density, pressure, and the x- and y-components of the magnetic field, are set to be outflow.

2.3. Implementing Ohmic and ambipolar diffusion

The temperature of the coronal plasma is typically of the order of 10⁶ K (see, e.g., Priest 2014), so the plasma is fully ionized. The temperature in prominences is much smaller, around 10⁴ K or lower (see, e.g., Tandberg-Hanssen 1995; Parenti 2014). As a consequence of that, the prominence plasma is only partially ionized (see, e.g., Labrosse et al. 2010). In a partially ionized plasma, ambipolar diffusion, caused by ion-neutral collisions, is an important dissipation mechanism (see, e.g., Osterbrock 1961; Pandey et al. 2008; Soler et al. 2009; Khomenko & Collados 2012; Ballester et al. 2018; Nóbrega-Siverio et al. 2020). Similarly, the Ohmic diffusion, caused by the collisions of electrons with other species, is heavily enhanced by the presence of neutrals. Both Ohmic and ambipolar diffusion can dissipate Alfvén waves in the partially ionized prominence medium (see, e.g., Soler et al. 2009).

In the presence of Ohmic and ambipolar diffusion, the resistive term $\widehat{\eta }·\mathit{j}$ in Eqs. (6) and (7) can be decomposed as (see, e.g., Bittencourt 2004),

$$\begin{array}{c}\hfill \widehat{\eta }·\mathit{j}=({\eta }_{\mathrm{O}}+{\eta }_{\mathrm{A}}{|\mathit{B}|}^2)\mathit{j}-{\eta }_{\mathrm{A}}(\mathit{B}·\mathit{j})\mathit{B},\end{array}$$(11)

where η_O and η_A are the Ohmic and ambipolar diffusion coefficients, respectively. The first term on the right-hand side of Eq. (11) only includes diagonal elements in the resistivity tensor, $\widehat{\eta }$. This term depends on both η_O and η_A. Conversely, the second term introduces both diagonal and off-diagonal elements and depends on η_A alone. In the PLUTO code, the resistivity tensor is assumed to be diagonal (Eq. (8)). The off-diagonal terms associated with the second term on the right-hand side of Eq. (11) cannot be included in the code. Therefore, the implementation of the ambipolar diffusion in PLUTO can only be done in an approximate manner (see, e.g., Khomenko & Collados 2012; Nóbrega-Siverio et al. 2020; Moreno-Insertis et al. 2022, for general implementations of the ambipolar diffusion). To circumvent this limitation of the code, we exploited the fact that, in the simulations, the background magnetic field along the tube axis is strong compared with the perturbations across the background magnetic field (see details in the Appendix A). Hence, in Eq. (11) we write

$$\begin{array}{c}\hfill \mathit{B}=B_0\widehat{z}+{\mathit{B}}_1,\end{array}$$(12)

where B₁ represents the perturbations over the background axial field and verify that |B₁|≪B₀. Consequently, we approximate $|\mathit{B}{|}^2=B_0^2$ and $(\mathit{B}·\mathit{j})\mathit{B}\approx B_0^2\mathit{j}·\widehat{z}$. Under this approximation, $\widehat{\eta }$ becomes a diagonal tensor as:

$$\begin{array}{c}\hfill \widehat{\eta }\approx \left(\begin{array}{ccc}{\eta }_{\mathrm{O}}+B_0^2{\eta }_{\mathrm{A}}& 0& 0\\ 0& {\eta }_{\mathrm{O}}+B_0^2{\eta }_{\mathrm{A}}& 0\\ 0& 0& {\eta }_{\mathrm{O}}\end{array}\right)=\left(\begin{array}{ccc}{\eta }_{\mathrm{C}}& 0& 0\\ 0& {\eta }_{\mathrm{C}}& 0\\ 0& 0& {\eta }_{\mathrm{O}}\end{array}\right),\end{array}$$(13)

where ${\eta }_{\text{C}}={\eta }_{\text{O}}+B_0^2{\eta }_{\text{A}}$ is the Cowling’s diffusion coefficient, which combines both the ambipolar and Ohmic diffusion coefficients. The Cowling diffusion accounts for the total magnetic diffusion across the magnetic field. We implemented in the PLUTO code this approximate resistivity tensor. Some tests about the numerical implementation can be checked in the Appendix B.

The Ohmic, η_O, and ambipolar, η_A, diffusion coefficients depend upon the plasma local properties and are calculated in PLUTO in the whole domain at each time step. The expressions to compute both coefficients are (see, e.g., Ballester et al. 2018):

$$\begin{array}{cc}& {\eta }_{\mathrm{O}}=\frac{{\alpha }_{\mathrm{e}}}{{\mu }_0{e}^2{n}_{\mathrm{e}}^2},\hfill \end{array}$$(14)

$$\begin{array}{cc}& {\eta }_{\mathrm{A}}=\frac{{\xi }_{\mathrm{n}}^2}{{\mu }_0{\alpha }_{\mathrm{n}}},\hfill \end{array}$$(15)

where ξ_n = 1 − ξ_i is the neutral fraction, e is the electron charge, and n_e is the electron number density. Moreover, α_e and α_n are the total friction coefficients for electrons and neutrals, respectively, which account for the collisions that these species have with all other particles. In our particular case, the species can be electrons, protons, or hydrogen atoms because we used a pure hydrogen plasma. The total friction coefficient of a species β can be obtained as the sum of the individual friction coefficients between different species, namely,

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

Following Braginskii (1965) and Spitzer (1968), the friction coefficient between two different charged species is:

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

where n_β (n_β′) is the number density of the species β (β′), m_ββ′ is the reduced mass defined as m_βm_β′/(m_β + m_β′), k_B is the Boltzmann’s constant, and ϵ₀ is the vacuum electrical permittivity. Moreover, lnΛ_ββ′ is the Coulomb logarithm, which is included to consider the ineffectiveness of the interactions between the charged species after some distance. According to Spitzer (1962; see also Vranjes & Krstic 2013), the Coulomb logarithm is:

$$\begin{array}{c}\hfill ln{\mathrm{\Lambda }}_{\beta \beta \prime }=ln\left[\frac{24\pi {\left({ϵ}_0{k}_{\mathrm{B}}T\right)}^{3/2}}{e^3\sqrt{n_{\beta }+n_{\beta \prime }}}\right]·\end{array}$$(18)

On the other hand, the friction coefficient between a charged and a neutral species is:

$$\begin{array}{c}\hfill {\alpha }_{\beta \mathrm{n}}=n_{\beta }{n}_{\mathrm{n}}\sqrt{\frac{8k_{\mathrm{B}}T}{\pi m_{\beta \mathrm{n}}}}{\sigma }_{\beta \mathrm{n}},\end{array}$$(19)

where the sub-index n denotes a neutral and σ_βn is the collisional cross-section, which we obtained from Vranjes & Krstic (2013). The values of the collisional cross-section are weakly dependent on the temperature (see Vranjes & Krstic 2013). However, for typical prominence temperatures considered here, this result can be safely ignored. Thus, we used constant values of the collisional cross-sections for simplicity.

We studied the values of the temperature, the ionization fraction, and the Ohmic, ambipolar, and Cowling diffusion coefficients in a longitudinal cut at the tube axis, r = 0, corresponding to the background model. The results are shown in Fig. 4. We found the minimum of temperature, T ≈ 8000 K, and the minimum of the ionization ratio, ξ_i ≈ 0.547, at the core of the prominence thread. The ambipolar diffusion is dominant for typical conditions of prominence threads, so it is the main contributor to the Cowling diffusion (see Melis et al. 2021, 2023). A similar situation occurs in the high chromosphere (Khomenko et al. 2014b). Magnetic diffusion is highly anisotropic in the partially ionized prominence plasma, being much more efficient across than along the magnetic field direction.

Fig. 4. Longitudinal cuts along the tube axis, r = 0, of the (a) Ohmic diffusion coefficient, (b) ambipolar diffusion coefficient, (c) Cowling diffusion coefficient, (d) temperature, and (e) ionization fraction in the equilibrium.

3. Numerical simulations

We aim to study the nonlinear evolution of standing torsional oscillations in the prominence thread. Unless otherwise stated, all quantities are specified in normalized units. Density, velocity, and length are normalized with respect to the external density, ρ_e, external Alfvén speed, v_A, e, and the thread radius, R, respectively, whose values have been given before. In turn, the normalized time, $\overline{t}$, is expressed in units of the Alfvén travel time, namely, t_A = R/v_A, e ≈ 1.53 s.

Before analyzing the result of the simulations, it is useful to calculate the expected period of the standing oscillations. Following the analytical treatment in Díaz-Suárez & Soler (2021b), to calculate the period of the longitudinally fundamental torsional Alfvén mode we solved the 1D wave equation for the linear azimuthal component of velocity, $v_{\phi }^{\prime }$, in the equilibrium, namely,

$$\begin{array}{c}\hfill \frac{{\partial }^2{v}_{\phi }\prime }{\partial t^2}=v_{\mathrm{A},0}^2(r,z)\frac{{\partial }^2{v}_{\phi }\prime }{\partial z^2}·\end{array}$$(20)

Equation (20) is transformed into an eigenvalue problem by expressing the temporal dependence as exp(iω_At) where ω_A is the Alfvén eigenfrequency. Thus,

$$\begin{array}{c}\hfill {\omega }_{\mathrm{A}}^2{v}_{\phi }\prime =-v_{\mathrm{A},0}^2(r,z)\frac{{\partial }^2{v}_{\phi }\prime }{\partial z^2},\end{array}$$(21)

with the boundary conditions of $v_{\phi }^{\prime }(r,z=\pm L/2)=0$. Because of the spatial dependence of v_A, 0, we solved the eigenvalue problem numerically with Wolfram Mathematica for the longitudinally fundamental mode. The period of the oscillations is P = 2π/ω_A. The period varies with the radial position because of the radial dependence of v_A, 0. The period in the uniform core of the thread is ∼535 t_A (13.62 min in physical time), while that in the corona is 100 t_A (2.54 min in physical time). This means that the Alfvén modes in the external plasma would have completed more than 5 periods of oscillation before the internal mode had completed a single period. Of course, the period continuously varies in the transition between the internal and external plasmas. The result of this continuum of periods will be a strong phase mixing in the radial direction.

3.1. Ideal MHD dynamics

First, we studied the dynamics of the oscillations in the absence of resistivity, so that we set η_O = η_A = 0 and performed an ideal MHD simulation. After the initial excitation, the flux tube is let to evolve. To help us understand better the dynamics of the simulation, besides studying the evolution of the density and the temperature, we also studied the vorticity, ω = ∇ × v, and the current density, $\mathit{j}={\mu }_0^{-1}\mathrm{\nabla }\times \mathit{B}$ because these quantities are extremely sensitive to spatial gradients in the velocity field and in the magnetic field, respectively. They can help us visualize the development of turbulence. The results of this ideal simulation are displayed in Figs. 5 and 6 in two different planes of the complete 3D box.

Fig. 5. Cross-sectional cuts of density (top left), temperature (top right), current density squared (bottom left), and vorticity squared (bottom right) at the end of the ideal MHD simulation, $\overline{t}=820$. Temperature is normalized with respect to the initial value at r = 0. Logarithmic scale is used for each variable except for the density to optimize visualization. The cross-sectional cuts are done at the tube center, z = 0, except for current density squared which is done at z = R. The complete temporal evolution is available as an online movie.

Fig. 6. Same as Fig. 5, but the cut is done longitudinally at y = 0. Note: the horizontal and vertical axes are not to scale. The complete temporal evolution is available as an online movie.

On the one hand, Fig. 5 shows cross-sectional cuts of the density (top left panel), the temperature (top right panel), the current density squared (bottom left panel), and the vorticity squared (bottom right panel). We used logarithmic scale to optimize the visualization for each variable except for the density. The cross-sectional cuts for density, vorticity squared, and temperature are done at the tube center, z = 0, while that for current density squared is done at z = R, which is a plane slightly displaced from the center. The reason for showing the current density in a different plane is that this variable has a node at the tube center. Since the relevant dynamics occurs in the core and in the transition layer of the prominence thread, we only showed a subdomain of the computational box where x, y ∈ [ − 1.5R, 1.5R]. We checked that nothing of interest happens outside this subdomain. The complete temporal evolution can be seen in the accompanying animation, while the still image in Fig. 5 displays the results at the final simulation time, $\overline{t}=820$, corresponding to ∼21 min in physical time.

On the other hand, Fig. 6 shows the same variables as in Fig. 5, but in a longitudinal plane to the tube corresponding to y = 0. As before, the results are only displayed in a subdomain of the whole box where x ∈ [ − 1.8R, 1.8R], while the z-direction is fully shown. In Fig. 6 we used a larger scale for the current density than in Fig. 5 to avoid early saturation and ease the visualization. Again, the complete temporal evolution can be seen in the accompanying animation, while the still image in Fig. 6 displays the results at the final simulation time.

As expected, we found that the torsional Alfvén modes oscillate with different frequencies in adjacent radial positions owing to the transversely nonuniform density. The oscillations become out of phase as time increases. Such a phenomenon is called phase mixing and is of linear nature (see, e.g., Heyvaerts & Priest 1983; Nocera et al. 1984; De Moortel et al. 2000; Smith et al. 2007; Prokopyszyn et al. 2019; Ebrahimi et al. 2020; Van Damme et al. 2020; Ebrahimi & Soler 2022). Due to phase mixing, the azimuthal component of velocity (not shown here) alternates between negative and positive values in adjacent radial positions. Such alternation generates azimuthal shear flows (see, e.g., Fig. 2 in Heyvaerts & Priest 1983 or Fig. 6 in Díaz-Suárez & Soler 2021b). These shear flows manifest themselves in the cross-sectional cuts of the vorticity and the current density as concentric rings, which are better seen in the animation of Fig. 5. Similar ring-shaped structures can also be found in simulations of torsional oscillations of coronal loops (Díaz-Suárez & Soler 2021b, 2022). These early structures in vorticity and current density are also equivalent to those seen in simulations of kink oscillations of coronal loops (Antolin et al. 2014; Howson et al. 2017b; Antolin & Van Doorsselaere 2019). However, the spatial distribution in the case of kink oscillations is not in the form of concentric rings owing to the different azimuthal symmetry that torsional and kink modes have. Besides, resonant absorption does not happen in our simulation, which is a concurrent process to phase-mixing in simulations of kink modes.

The annular structures in vorticity and current density caused by phase mixing are seen in the simulations before any appreciable disturbance in the density or the temperature. At the early stages of the evolution, we can see some weak perturbations in density in the form of periodic accretion (voiding) of plasma to (from) the center of the thread, which is better appreciated in the animation of Fig. 6. These density variations are not related to phase mixing but are caused by the ponderomotive force, a nonlinear effect that couples Alfvén and slow magnetosonic modes (Hollweg 1971; Rankin et al. 1994; Tikhonchuk et al. 1995; Terradas & Ofman 2004).

As phase mixing develops, the azimuthal shear flows gradually intensify and eventually trigger the KHI (see, e.g., Heyvaerts & Priest 1983; Browning & Priest 1984; Soler et al. 2010; Zaqarashvili et al. 2015; Barbulescu et al. 2019). As these flows are perpendicular to the background magnetic field lines, the magnetic tension cannot avoid the triggering of KHI (see, e.g., Chandrasekhar 1961). We can visually identify in Fig. 5 that the first KHI-associated deformations of temperature and density occur at τ_vis ≈ 280 (∼7 min in physical time), which corresponds to less than three periods of the external torsional Alfvén mode and about half a period of the internal mode. However, the examination of the form of the vorticity for that time already reveals the existence of significant deformations in the boundary between the inhomogenous layer and the external medium. This makes us wonder whether the visual estimation of the KHI onset time from the density evolution might overestimate the actual onset time. Indeed, the initial growth of KHI perturbations can be appreciated in the vorticity at an earlier time than in the density and the temperature. A more detailed analysis of onset time is given in Sect. 3.1.1.

Once the KHI is triggered, large eddies are formed and rapidly grow inside the nonuniform boundary layer of the tube. The KHI evolves nonlinearly, breaking the initially large eddies into small eddies leading naturally to turbulence. Turbulence mixes the hot and light coronal plasma with the cold and dense prominence thread plasma. Initially, turbulence occurs locally in the nonuniform boundary layer alone, but as time increases, part of the core and the external medium are also affected. In the temperature evolution, we can see intrusions of hot plasma towards the cool core of the thread and, in density, an apparent breakup of the prominence material. Importantly, neither the KHI nor turbulence are seen in the longitudinal direction (see Fig. 6), meaning that the turbulence only develops perpendicularly to the magnetic field lines. In the case of torsional oscillations of coronal loops, this matter has been investigated in Díaz-Suárez & Soler (2021b); Díaz-Suárez & Soler (2022).

3.1.1. Estimating the Kelvin–Helmholtz instability onset time

To estimate quantitatively the KHI onset time, we used the discrete Fourier transform in the azimuthal direction of the azimuthal and radial velocity components. The KHI is expected to excite high azimuthal modes. This technique was first used by Terradas et al. (2018) and later by Antolin & Van Doorsselaere (2019) and Díaz-Suárez & Soler (2021b, 2022) in coronal loops. Particularly, we investigated the azimuthal and radial components of velocity in a cross sectional cut at the tube center, z = 0, and in the boundary between the external medium and the inhomogeneous density layer, that is, for r = 1.3R. This radial position is chosen because we visually identify the growth of the first KHI eddies at that location. We applied the discrete Fourier transform to both velocity components using the fast Fourier transform (FFT) algorithm with the Scipy module (Virtanen et al. 2020). Following the same notation as in Terradas et al. (2018), the discrete Fourier transform can be set as:

$$\begin{array}{c}\hfill G(p)={\displaystyle \sum _{k=0}^{N-1}}g(k)exp(-\frac{2\pi ipk}{N})·\end{array}$$(22)

In Eq. (22), N is the number of samples, g(k) is the angular sampling of the azimuthal (radial) velocity, and p is an integer that plays the role of the azimuthal wavenumber and ranges between 0 and N − 1. In this notation, p = 0 is the torsional or sausage mode, p = 1 is the kink mode, and p ≥ 2 are the fluting modes (see, e.g., Roberts 2019). Generally, the Fourier amplitudes of Eq. (22) are complex with the exception of G(p = 0). Consequently, we calculated the modulus of all the Fourier coefficients. The temporal evolution of the Fourier amplitudes is shown in Fig. 7. To optimize the visualization of the contribution of the modes different from p = 0, which would be the dominant one, we added the first twenty modes starting from p = 1.

Fig. 7. Top panel: temporal evolution of the Fourier coefficient with p = 0 (red solid line) and of the sum of the modulus of the first twenty positive Fourier coefficients (blue dashed line). For this analysis, the azimuthal component of velocity along a circle of radius 1.3R at the z = 0 plane is used. The two vertical dashed black lines show, respectively, the onset of the KHI obtained from this Fourier analysis, τKH ≈ 70, and the time at which the KHI is first seen to grow in the density, τvis = 280. Bottom panel: same as the top panel but for the radial component of velocity.

As expected, the p = 0 mode is dominant during most of the simulation because this is the azimuthal symmetry imposed by the initial excitation. The sum of the first twenty Fourier modes other than p = 0 is initially zero, as expected. The amplitude of the p = 0 mode oscillates with a periodicity of 100 in code units, which consistently matches the period of the local torsional Alfvén mode at the chosen radial position. As the linear development of phase mixing does not excite higher azimuthal modes, the time at which |G(p > 0)| departs from zero can be defined as the KHI onset time, τ_KH. In particular, we find that τ_KH ≈ 70 in code units. We recall that the periods of the internal and external torsional modes are 535 and 100 in code units, respectively, which indicate that the KHI is very quickly driven in the system. The KHI onset happens significantly earlier than the visually determined time of τ_vis ≈ 280 that is based on the growth of the KHI vortices in density. In fact, the sum of the first twenty Fourier modes is already comparable with the amplitude of the torsional mode when the KHI eddies become visible in density and in temperature.

A similar analysis but with the inclusion of higher azimuthal modes does not change the obtained results (not shown here). However, when the analysis is done in layers of the transition region nearer the core, the period of the torsional mode increases, which is expected, and the growth of the other Fourier modes is delayed. The chosen radial position of r = 1.3R corresponds to the location where the KHI first grow, so that the analysis at that location provides the smallest KHI onset time.

3.1.2. Evolution of integrated vorticity and current density

In the later stages of the simulation when turbulence plays a predominant role, we notice the very small scales that are generated in the current density and the vorticity. In addition, the values of the two quantities seem to increase with time during both the initial quasi-linear stage governed by phase mixing and the later nonlinear stage governed by turbulence. To quantify the increase, we calculated the vorticity squared and the current density squared integrated in the whole computational domain as functions of time, namely

$$\begin{array}{cc}& {\mathrm{\Omega }}^2(t)=\int {|\mathit{\omega }(\mathit{r},t)|}^2{\mathrm{d}}^3\mathit{r},\hfill \end{array}$$(23)

$$\begin{array}{cc}& I^2(t)={\displaystyle \int {|\mathit{j}(\mathit{r},t)|}^2}{\mathrm{d}}^3\mathit{r}.\hfill \end{array}$$(24)

The calculations are shown in Fig. 8 after normalizing Ω² with respect to the initial value. The background model is current-free, so that I² = 0 at t = 0. Using the vorticity evolution in the top panel of Fig. 8, we can distinguish three stages in our simulation: a quasi-linear phase, a nonlinear phase dominated by the KHI growth, and a saturation phase where turbulence develops. The two vertical black dashed lines in Fig. 8 denote the transitions between these distinct stages.

Fig. 8. Temporal evolution of the vorticity squared (top panel) and the current density squared (bottom panel) integrated in the whole computational domain. The vorticity values are normalized with respect to the initial value. The two dashed black lines show respectively the onset of the Kelvin–Helmholtz instability inferred from the analysis of azimuthal modes, τKH = 70, and the time at which the integrated vorticity squared saturates, τsat = 300.

In the short-lived quasi-linear phase, there is an approximately linear increase of Ω² with time owing to the slow building up of small scales by phase mixing (Soler & Terradas 2015; Howson et al. 2020; Díaz-Suárez & Soler 2021b). This stage ends when the KHI is triggered for t = τ_KH, which coincides with a change of trend in the Ω² curve. This coincidence reinforces the validity of the Fourier analysis to determine the KHI onset time. During the second phase, the KHI development more rapidly increases the values of vorticity (see Howson et al. 2017b; Guo et al. 2019; Díaz-Suárez & Soler 2021b, 2022). This large increase in vorticity should continue indefinitely in our ideal simulation as the large vortices break into smaller vortices during the KHI turbulent evolution. Nevertheless, the increase in Ω² stops and saturates at t = τ_sat, with τ_sat ≈ 300 in code units for this particular simulation. The saturated phase is characterized by the development of turbulence in the flux tube. The limited spatial resolution and its associated numerical diffusion is behind this saturation. Having a high spatial resolution is crucial to capture the finest structures generated by the turbulence (Howson et al. 2017a; Díaz-Suárez & Soler 2021b).

There is also an increase in the values of I², but much smoother than that observed for Ω². The three separate phases discussed before are not so clearly distinguishable in the evolution of I². One may ask why both variables have different behaviors. To answer that question, one needs to consider the spatial dependence of vorticity and current density along the flux tube. For a standing torsional Alfvén wave, the velocity perturbations have nodes at the tube footpoints and an antinode at the tube apex. The vorticity follows the same dependence. Therefore, in Eq. (23) the largest contribution to the integral of Ω² comes from the apex of the tube. In turn, the magnetic field perturbations have antinodes at the tube footpoints and a node at the tube apex, with the current density mimicking the same behavior. Hence, in Eq. (24) the largest contribution to the integral of I² originate from the footpoints of the magnetic flux tube. Since the KHI and turbulence develop predominantly around the tube apex, they heavily impact on the evolution of Ω². Conversely, the KHI and turbulence do not have such a pronounced development at the tube footpoints, where phase-mixing remains as the dominant mechanism that creates small scales over time. As a result, there is a less significant imprint of the KHI and turbulence in the evolution of I², which keeps displaying the linear growth caused by phase mixing even after Ω² has already saturated.

3.2. Effect of Ohmic and ambipolar diffusion

The generation of small structures in the current density and the increase of its value over time might lead to the dissipation of the wave energy if magnetic diffusion is included. Here we focus on studying the role of Ohmic and ambipolar diffusion. To this end, we ran two additional simulations. These new simulations were performed under the same conditions as the ideal simulation but adding the resistive term in the equations. In the first simulation, called the Ohmic simulation, we included Ohmic diffusion, but not ambipolar diffusion. Consequently, no approximation is used in the treatment of resistivity and the PLUTO resistivity tensor is in this case isotropic, namely

$$\begin{array}{c}\hfill \widehat{\eta }=\left(\begin{array}{ccc}{\eta }_{\mathrm{O}}& 0& 0\\ 0& {\eta }_{\mathrm{O}}& 0\\ 0& 0& {\eta }_{\mathrm{O}}\end{array}\right),\end{array}$$(25)

In the second simulation, called the Cowling simulation, we included both ambipolar diffusion and Ohmic diffusion, so the resistivity tensor is anisotropic and follows Eq. (13). The purpose in this subsection is to compare the two dissipative simulations with the previously discussed ideal MHD simulation.

The overall, large-scale dynamics of both dissipative simulations are identical to that of the ideal MHD simulation discussed in Sect. 3.1. No noticeable differences are seen in the results of density, temperature, and vorticity. We remind the reader that the values of η_O and η_A used in the dissipative simulations are the realistic ones in the partially ionized prominence plasma (see Fig. 4). As a consequence of that, dissipation is not artificially enhanced in our simulations and its influence is necessarily reduced to the smallest scales that develop during the evolution. The current density is the variable that displays more significant differences between simulations, but even for that variable, the differences can be considered as small. Figure 9 shows a cross-sectional at z = R of the current density squared in logarithmic scale for the ideal, Ohmic, and Cowling simulations. We find slight differences in the fine scales of the turbulence that develops at large times in the evolution, but no clear relation can be deduced by comparing the results of the three simulations. It is obvious that the presence of dissipation affects somehow the development of the smallest scales. This fact, together with the intrinsic chaotic nature of turbulence (see, e.g., Biskamp 2003) causes the appearance of slightly different turbulent patterns for the current density in the three simulations.

Fig. 9. Cross-sectional cut at z = R of the current density squared in logarithmic scale at the end of the simulations, $\overline{t}=820$. The panels are sorted as follows: ideal simulation (left), Ohmic simulation (center), and Cowling simulation (right). The complete temporal evolution is available as an online movie.

To quantify the differences in the current density between the three simulations, Fig. 10 shows the temporal evolution of the current density squared integrated in the whole computational domain. The curve corresponding to the ideal simulation was already displayed in the bottom panel of Fig. 8. The three cases behave similarly until turbulence begins to dominate the dynamics late in the evolution. The three curves separate from that time onwards, although an increasing trend remains in the three cases. This increasing trend points out that the generation of small scales in current density continues even when magnetic diffusion works to dissipate those small scales. Counter-intuitively, the ideal simulation shows lower values of current density than the two dissipative simulations in some time range. In principle, one should expect Ohmic and ambipolar diffusion to slow down the KHI development and inhibit the formation of the smaller KHI vortices compared to the ideal case. In this line of thought, the increase in the current density in the dissipative simulations should be less pronounced than in the ideal simulation. However, the results do not point in that direction. Instead, we find no clear pattern relating the efficiency of physical dissipation with the behavior of the integrated current density. The reason why physical dissipation is acting inefficiently to damp the perturbations in the current density resides in the realistically small values of the diffusion coefficients and the likely influence of the unavoidable numerical dissipation. Since we are not able to disentangle the role of numerical dissipation from that of physical dissipation, it is not possible for us to determine the relative importance of numerical dissipation. However, the fact that the three simulations show a different turbulent pattern at the very small scales indicates that numerical dissipation does not entirely dominates and that physical dissipation is playing a role.

Fig. 10. Same as bottom panel from Fig. 8 for the ideal simulation (solid red line), but now including the results from the Ohmic simulation (blue dashed line) and the Cowling simulation (black dash-dotted line).

We go on to examine the heating produced in the dissipative simulations. Using the approximate resistivity tensor of Eq. (13), the volumetric heating rate, Q, is:

$$\begin{array}{c}\hfill Q={\mu }_0(\widehat{\eta }·\mathit{j})·\mathit{j}\approx {\mu }_0(\eta |{\mathit{j}}_{\Vert }{|}^2+{\eta }_{\mathrm{C}}{|{\mathit{j}}_{\perp }|}^2),\end{array}$$(26)

where j_∥ and j_⊥ are the components of the current density in the parallel and perpendicular directions to the background magnetic field, calculated as:

$$\begin{array}{cc}& {\mathit{j}}_{\Vert }=\frac{\mathit{j}·\mathit{B}}{{|\mathit{B}|}^2}\mathit{B}\approx j_z,\hfill \end{array}$$(27)

$$\begin{array}{cc}& {\mathit{j}}_{\perp }=\frac{\mathit{B}\times (\mathit{j}\times \mathit{B})}{{|\mathit{B}|}^2}\approx j_x\widehat{x}+j_y\widehat{y}.\hfill \end{array}$$(28)

When ambipolar diffusion is absent, η_C = η, and the heating is isotropic. In the presence of ambipolar diffusion, the heating is expected to be more intense and dominated by the dissipation of perpendicular currents. Figure 11 shows the volumetric heating rate in the Ohmic and Cowling simulations at the same cross-sectional cut as that in Fig. 10.

Fig. 11. Snapshot of the volumetric heating rate, Q, in a cross sectional cut at z = R at the end of the simulation time, $\overline{t}=820$. A logarithmic scale has been used in both panels for an optimal visualization. Left panel: Ohmic simulation. Right panel: Cowling simulation. The complete temporal evolution is available as an online movie.

In the Ohmic simulation, we find that the heating is negligible in the initial stage of the dynamics governed by phase mixing alone. Once the KHI evolves nonlinearly and turbulence develops, there is a large increase in the values of the heating rate. The heating is essentially produced in the nonuniform transitional layer between the core of the prominence thread and the corona, where turbulence mainly develops. Negligible Ohmic heating is obtained in the cool core of the thread, since turbulence does not completely reach that region of the flux tube.

The heating obtained in the Cowling simulation is around three orders of magnitude larger than the heating in the Ohmic simulation. Such a result is consistent with the larger efficiency of ambipolar diffusion in prominence threads (see, e.g., Melis et al. 2021, 2023). Again, the heating is most important in the nonuniform transitional layer, but now a non-negligible heating is obtained in the cool core of the thread. The fact that the cool core is only partially ionized results in the presence of heating although turbulence does not develop in that region. Thus, the heating in the cool core is unrelated to turbulence and is directly caused by the ambipolar dissipation of the Alfvén waves (Soler et al. 2009).

The regions where heating is more intense appear to be very spatially localized in the turbulent annulus, which might limit the global effect that this heating could have on the prominence thread. To explore whether the obtained local heating has a sizeable impact on the thermodynamic state of the thread, we selected a subdomain of the model that encompasses the partially ionized part of the thread. The considered subdomain corresponds to x/R ∈ [−1.4,1.4], y/R ∈ [−1.4,1.4], and z/R ∈ [−4,4]. We studied the evolution in time of the average temperature in that region. The calculation is shown in Fig. 12. For comparison purposes, we also included the result corresponding to the ideal MHD simulation, for which there is no physical heating. Despite the fact that heating is present in the Ohmic and Cowling simulations, the average temperature in the three cases shows a similar behavior. During the initial, quasi-linear phase the three curves superimpose. In that phase, we find that the average temperature remains roughly constant, with some slight increases and decreases probably caused by the adiabatic compression and expansion of the plasma along the flux tube due to the ponderomotive force. When the KHI and the turbulence occur, the average temperature starts to decrease. The evolution of the average temperature in the dissipative simulations is practically the same as in the ideal simulation in this decreasing phase too. The average temperature decreases in all cases due to the mixing of the hot coronal plasma and the cold prominence plasma owing to the nonlinear evolution of the KHI and the turbulence, as Hillier & Arregui (2019) and Hillier et al. (2023) explained. However, we note that unlike in Hillier et al. (2023), in our simulations we do not include radiative losses, which would further cool the plasma until achieving a radiative equilibrium. The decrease of the average temperature that we find in our simulations is exclusively due to the mixing of the internal and external plasmas. At the end of the simulations, the average temperature in the considered subdomain has decreased about 2% with respect to the initial value, so that this apparent cooling of the thread due to the plasma mixing is indeed quite modest, but still much more important that the Ohmic and ambipolar heating. We found a similar evolution for the average temperature using larger subdomains in z- x- and y-directions, as long as the mixing layer is included (not shown here).

Fig. 12. Temporal evolution of the average temperature in a subdomain where x/R ∈ [−1.4,1.4], y/R ∈ [−1.4,1.4], and z/R ∈ [−4,4] for the ideal simulation (solid red line), the Ohmic simulation (blue dashed line) and the Cowling simulation (black dash-dotted line). The temperature is normalized with respect to the initial average temperature.

We conclude that the plasma heating in the dissipative simulations does not play a relevant role in the thermodynamic state of the thread. The heating is very localized in a turbulent annulus around the cool core of the thread and its contribution is completely overwhelmed by the effective cooling caused by the plasma mixing. The results here are in apparent contradiction to those of Melis et al. (2021), who found that heating by Alfvén waves in prominence threads could be important enough to balance energy losses due to radiative cooling. However, there are important differences between the work of Melis et al. (2021) and the present work, so that their results cannot directly be compared with our findings. Melis et al. (2021) studied propagating waves in a frequency range far higher than in our work (they considered frequencies as high as 1 Hz), while here we excited standing waves that have much lower frequencies (between 1.22 mHz and 6.54 mHz). The efficiency of Alfvén wave dissipation increases with the frequency, so that the propagating waves studied by Melis et al. (2021) are more efficiently damped than the standing modes studied here. In addition, Melis et al. (2021) used a 1D thread model and neither the KHI nor the turbulence are able to develop in their configuration.

4. Forward modeling

Here we explore the possible observational signatures that the nonlinear evolution of torsional Alfvén oscillations in prominence threads may leave. To this end, and inspired by the work of Martínez-Gómez et al. (2022) in the case of kink oscillations, we computed synthetic Hα observations,

We used the method of fast synthesis of the Hα line profile given in Heinzel et al. (2015). The details of the method can be found in Heinzel et al. (2015) and we only give a summary of its implementation. The formal solution to the radiative transfer equation in the line-of-sight (LOS) direction is:

$$\begin{array}{c}\hfill I_{\nu }=S[1-exp(-{\tau }_{\nu })],\end{array}$$(29)

where S is the source function assumed uniform and τ_ν is the optical thickness in the LOS, defined as:

$$\begin{array}{c}\hfill {\tau }_{\nu }={\displaystyle {\int }_{\mathrm{LOS}}}\kappa (\nu ,l)\mathrm{d}l,\end{array}$$(30)

where l is the LOS coordinate and κ(ν, l) is the absorption coefficient given by

$$\begin{array}{c}\hfill \kappa (\nu ,l)=\frac{\pi e^2}{m_{\mathrm{e}}c}{f}_{23}{n}_2(l)\varphi (\nu ,l).\end{array}$$(31)

In Eq. (31), f₂₃ is the Hα line oscillator strength (see, Goldwire 1968, for the tabulated value), c is the speed of light, and m_e is the electron mass. Moreover, n₂ is the population of hydrogen atoms in the second quantum level and ϕ(ν, l) is the normalized absorption profile, which is approximated by a Gaussian profile (see Heinzel et al. 2015). From here on, we assume that the LOS direction is the y-direction, so the x- and z-directions form the plane of sky (POS). Bearing in mind the local Doppler shifts owing to the LOS velocity in the y-direction, we write,

$$\begin{array}{c}\hfill \varphi (\nu ,y)=\frac{1}{\sqrt{\pi }\mathrm{\Delta }{\nu }_{\mathrm{D}}(y)}exp\{-\frac{{[\nu -{\nu }_0{(1+v_y(y)/c)}^{-1}]}^2}{\mathrm{\Delta }{\nu }_{\mathrm{D}}^2(y)}\},\end{array}$$(32)

where ν₀ is the Hα rest frequency, v_y(y) is the y-component of velocity, and Δν_D(y) is the thermal and microturbulent broadening given in Heinzel et al. (2015). The microturbulent broadening depends on the microturbulent velocity, whose value is set to 5 km s⁻¹ as in Heinzel et al. (2015). In the simulations, we occasionally found maximum LOS velocities of ±10 km s⁻¹, which is the limit of accuracy of the method of Heinzel et al. (2015). However, the typical LOS velocities are normally around ±6 km s⁻¹, which are compatible with observationally reported values (see, e.g., Schmieder et al. 2010; Gunár et al. 2012).

The hydrogen second-level population, n₂, can be calculated from the electron density. Following Heinzel et al. (2015), both physical quantities are related as

$$\begin{array}{c}\hfill n_2(y)=\frac{n_{\mathrm{e}}^2(y)}{f(T(y),p(y))},\end{array}$$(33)

where f(T(y),p(y)) is a function that depends on temperature and gas pressure. Physically, the function f is associated with the rate of the photoionization and the radiative recombination from and to the hydrogen second-level population (Heinzel et al. 1994). We calculated the values of the function f using bilinear interpolation from the values in Table 1 in Heinzel et al. (2015), assuming a constant altitude of 10 Mm. We are aware that the source function can vary with height (see, e.g., Heinzel et al. 1994). However, for consistency with our choice, we neglected the variation of the source function with height.

Since the ideal and dissipative simulations are nearly identical in the description of the overall dynamics of the thread, we only use here the results from the ideal simulation to calculate the synthetic observation. The specific intensity at the Hα rest frequency in the POS is displayed in Fig. 13 and its accompanying animation. As the fully ionized plasma does not emit in Hα, we used a reduced numerical domain in z where z/R ∈ [ − 10, 10]. The emission occurs at the core of the partially ionized prominence thread, which is the only part of the model that is visible in the synthetic images. In the temporal evolution, first one can see a periodic pulsation in the intensity that emerges from the thread. This initial phase is then followed by the appearance of fine strands that are longitudinal to the thread axis and have a brighter intensity than the rest of the thread. The reason for the existence of these two distinct phases in the synthetic observation is analyzed below.

Fig. 13. Normalized specific intensity at the Hα rest frequency in four different times. The normalization is done with respect to the maximum value at the beginning of the simulation. We note that the longitudinal, z, and transverse, x, directions to the thread are not to scale. A complete temporal evolution is available as an online movie.

The initial periodic variations in the Hα intensity are caused by the integration of the LOS flows and, to a lesser extent, the ponderomotive flows along the thread. To confirm this hypothesis, first, we created a time-distance map of the Hα intensity at z = 0, that is, across the thread. This is displayed in Fig. 14 (left panel). The initial intensity pulsations are clearly visible. The period of such pulsations is ≈6.6 min and roughly two periods are seen in the time-distance map. Then, we repeated the time-distance map but neglecting the LOS flows. This is displayed in Fig. 14 (right panel). The comparison of the two maps reveals that the LOS flows or, more precisely, the integration of those flows along the LOS, have the net effect of producing the intensity pulsations. In the initial stages of the evolution, the LOS flows are caused by the phase mixing of the Alfvén modes. Therefore, the intensity pulsations appear to be an observational manifestation of phase mixing. To the best of our knowledge, this has not been discussed before in the literature.

Fig. 14. Time-distance map of the normalized specific intensity across the thread at the rest frequency of Hα (left). The normalization is done with respect to the maximum value at the beginning of the simulation. Right panel is the same as the left, but neglecting the LOS velocities.

Simultaneously to the pulsations, one can see a gradual brightening at the core of the prominence thread. This brightening is present in both panels of Fig. 14, although it is visually more evident in the right panel where the pulsations are absent. This gradual brightening is caused by the ponderomotive force, which tends to accumulate plasma around the thread center so increasing its density. Using Eq. (14) from Díaz-Suárez & Soler (2021b), we obtained the period of the associated slow MHD wave due to the ponderomotive force to be ≈24.5 min. To compute this period, we considered the longitudinally averaged density along the axis of the flux tube. The ponderomotive force acts so slowly that the KHI develops before a complete period of the slow MHD wave. In the time-distance map, the KHI and the onset of turbulence are seen in the form of an irregular pattern of bright patches ultimately caused by the development of the KHI vortices and the mixing of plasma. We note that the omission of the LOS flows does not alter much the structure of the bright patches, which confirms that their origin resides in the turbulent distribution of the density.

By the time that the KHI develops after ∼12 min, Fig. 13 shows the formation of fine strands at the core of the prominence thread. These strands are the observational signature of the KHI vortices, which develop perpendicularly to the tube axis. This is not a new result. The formation of strands owing to the KHI has been found in numerical simulations of magnetic tubes oscillating in the kink mode (e.g., Antolin et al. 2014, 2015, 2016; Antolin et al. 2017; Guo et al. 2019; Shi et al. 2021; Martínez-Gómez et al. 2022). Most of the previous works performed forward modeling of coronal loops, which involves spectral lines associated to much hotter temperatures than those of prominence threads.

Martínez-Gómez et al. (2022) studied the Hα emission of prominence threads oscillating in the kink mode and our results can be compared to theirs. Unlike in Martínez-Gómez et al. (2022), we do not see a global lateral oscillation of the flux tube, which is characteristic of kink modes. Another difference is that the initial periodic pulsations inside the prominence thread due to the phase-mixing flows along the LOS are not reported in Martínez-Gómez et al. (2022). Regarding the formation of the strands, Martínez-Gómez et al. (2022) found that the width of the strands in their kink mode simulations is between 10 km and 125 km. We found that the width of the strands in our torsional mode simulations range between 115 km and 168 km, approximately. Considering such size of the strands, they might be observable with instruments such as CRisp Imaging SpectroPolarimeter (CRISP; Scharmer et al. 2008) at the Swedish 1 m Solar Telescope (SST; Scharmer et al. 2003), or the Visible Image Spectrometer (VIS; Cao et al. 2010) at the Goode Solar Telescope (GST; Goode et al. 2010). Both instruments has an spatial resolution approximately equal to 0.1″ (∼70 km).

However, we should note that the radius of our simulated thread, 1 Mm, is larger than the typically reported thread radii, ∼200 − 300 km (see, e.g., Lin et al. 2005, 2007). An average value of observed thread radii is 228 km (see, e.g., Arregui et al. 2018), which is a factor of 4.38 smaller than the radius of the simulated thread. After scaling the width of the strands by the same factor, we obtain that the strands should have an actual width between 26 km and 38 km, approximately. Such fine strands cannot be seen with current instruments. Conversely, the next generation of solar telescopes might be able to observe these fine strands. The Visible Broadband Imager (VBI; Wöger et al. 2021) and the Visible Tunable Filter (VTF; Schmidt et al. 2014), both installed at Daniel K. Inouye Solar Telescope (DKIST; Rimmele et al. 2020), have spatial resolutions of 20 km and ∼25 km, respectively. The Tunable Imaging Spectropolarimeters (TIS) installed at the European Solar Telescope (EST; Quintero Noda et al. 2022) have an spatial resolution of ∼29 km at the rest wavelength of Hα in air. We repeated the synthetic modeling considering different LOS directions, but all cases showed similar results.

5. Concluding remarks

In this numerical work, we study the nonlinear evolution of standing torsional Alfvén waves in a low-β prominence-thread model embedded in a constant and axial magnetic field. The model consists of a flux tube that has a core, a transversely nonuniform transition region, and a external medium. The external medium has a uniform temperature and density while the longitudinally nonuniform core of the prominence thread is denser and cooler than the external medium. Both regions are continuously connected everywhere through a nonuniform transition. Moreover, the magnetic field, where the prominence thread is embedded, is line-tied at the photosphere. For simplicity, no chromospheric layer is included.

We excited the longitudinally fundamental mode of standing torsional Alfvén waves and performed three simulations using the PLUTO code: an ideal MHD simulation and two nonideal simulations including Ohmic diffusion alone and Ohmic together with ambipolar diffusion. Other nonideal effects, such as the Hall effect and the Biermann battery effect, were ignored as they are of much less relevance in prominences (Khomenko et al. 2014a; Ballester et al. 2018; Melis et al. 2021). The simulated dynamics undergoes three differentiated stages. The first stage is quasi-linear and dominated by the phase mixing of the Alfvén modes (see, e.g., Heyvaerts & Priest 1983). The second stage begins when phase-mixing azimuthal flows trigger the KHI. Since the KHI excites higher azimuthal modes, the KHI onset time has been obtained through an analysis of Fourier modes in the azimuthal direction using the azimuthal and radial components of the velocity at a specific radial position. The onset time obtained from the Fourier analysis agrees with the time at which a change of trend in the integrated values of vorticity happens. Once the KHI is triggered, large eddies initially appear, which evolve nonlinearly breaking into smaller and smaller eddies. This dynamics eventually leads to turbulence, so the dynamics of the thread is similar to that obtained in simulations of torsional oscillations of coronal loops (see Díaz-Suárez & Soler 2021b). Nonetheless, the limited spatial resolution of our ideal MHD simulation imposed a third phase, namely the saturated phase. In this stage, the integrated vorticity remains roughly constant when it should increase indefinitely. This situation, however, does not impede that the generation of turbulence is captured by our simulations.

Overall, the large-scale dynamics in the two dissipative simulations is nearly identical to that of the ideal simulation. Some differences appear only at the very small scales, specially in the current density. This rather subtle effect of dissipation, indeed, causes some impact in the evolution of turbulence, but we find it difficult to disentangle the effect of the physical dissipation from that of the inherent numerical dissipation. The plasma heating associated with the dissipation of currents is weak and very localized in an annulus region at the thread boundary. As a result of that, no global heating of the prominence thread is produced. Instead, the mixing of the cool prominence plasma with the hot coronal plasma leads to a moderate cooling of the thread towards an intermediate temperature (see Hillier et al. 2023). Thus, Ohmic and ambipolar dissipation of turbulence induced by torsional waves does not appear to be a mechanism capable of heating prominence threads.

We wondered what observational signatures the KHI and the turbulence might leave. To this end, we followed the method of Heinzel et al. (2015) to compute synthetic Hα observations. Before the onset of the KHI, we found a periodic brightening of the Hα intensity at the core of the prominence thread. This periodic brightening is caused by the integration of the flows along the LOS and, to a lesser extent, the longitudinal flows driven by the ponderomotive force. Later, we saw the formation of fine bright strands parallel to the thread axis owing to the KHI vortices and the turbulence. Forward modeling of coronal loop kink oscillations (see, e.g., Antolin et al. 2014, 2016, 2017; Guo et al. 2019; Shi et al. 2021) and prominence thread kink oscillations (Martínez-Gómez et al. 2022) reported similar strand formation. We showed that standing torsional oscillations are also able to drive such fine strand formation. When scaled according to the considered radius of the thread in the simulations, the effective width of the strands ranges between 26 km and 38 km, approximately. The next generation of solar telescopes DKIST (Rimmele et al. 2020) and EST (Quintero Noda et al. 2022) might be able to observe them. We note that while we modeled a prominence thread as a straight magnetic tube, in reality the prominence material is deposited in dips of the magnetic field due to the effect of gravity. Thus, their inclusion in future works would be more realistic. The study of the interaction between neighboring threads in a prominence would also be an interesting extension of this work.

Movies

Movie 1 associated with Fig. 5 (fig05_multipanel) Access here

Movie 2 associated with Fig. 6 (fig06_longcuts_def) Access here

Movie 3 associated with Fig. 9 (fig09_jsquaredatz1) Access here

Movie 4 associated with Fig. 11 (fig11_logheatz1def) Access here

Movie 5 associated with Fig. 13 (fig13_halpha) Access here

Acknowledgments

This publication is part of the R+D+i project PID2020-112791GB-I00, financed by MCIN/AEI/10.13039/501100011033. S.D.S. acknowledges the financial support from MCIN/AEI/10.13039/501100011033 and European Social Funds for the predoctoral FPI fellowship PRE2018-084223. This research was supported by the International Space Science Institute (ISSI) in Bern, through ISSI International Team project #457 (The Role of Partial Ionization in the Formation, Dynamics and Stability of Solar Prominences). We acknowledge the UIB for the use of the Foner cluster. S.D.S. thanks L. Melis, M. Kriginsky, and F. Bailén Martínez for their help. We also thank the anonymous referee for the constructive comments that improved the quality of the paper. During the analysis of data, we have used VisIT (Childs et al. 2012), Mathematica (Wolfram Research Inc. 2019), and Python 3.6. The Python packages that we have used are Matplotlib (Hunter 2007), Scipy (Virtanen et al. 2020) and Numpy (Harris et al. 2020). We are thankful to B. Vaidya and his contributors for the tool pyPLUTO.

References

Appendix A: Justification for neglecting off-diagonal terms in the resistivity tensor

As mentioned in Sect. 2.3, the resistivity tensor is defined in PLUTO as a diagonal tensor (see Eq. (8)). However, the ambipolar diffusion introduces off-diagonal elements, as shown in Eq. (11). Nonetheless, the properties of our problem allows us to implement the ambipolar diffusion in an approximate manner. Particularly, the background magnetic field, that is aligned with the main axis of the flux tube along the z-direction, is much stronger than the components across the background magnetic field, namely B_x and B_y. This observation allows us to neglect the off-diagonal elements associated with ambipolar diffusion. For the ideal MHD simulation, we computed the maximum values of B_x/B_z and B_y/B_z and plotted them against time. These results are displayed in Fig. A.1. As one can see, the maximum values of the magnetic field components perpendicular to the background field are less than 5% the value of the background field strength. In addition, we checked that these maximum values typically occur near the footpoints of the magnetic tube, where the plasma is fully ionized and ambipolar diffusion is absent. Therefore, the approximation explained in Sect. 2.3 to implement the ambipolar diffusion in the particular case under study is justified.

Fig. A.1. Temporal evolution of the maximum value of Bx/Bz (top). The same, but for By/Bz (bottom). The results are obtained from the ideal MHD simulation.

Appendix B: Test implementation of Ohmic and ambipolar diffusion in PLUTO

B.1. For 1.5 dimensions

Following Balsara (1996), Soler et al. (2009), or Ballester et al. (2018), the dispersion relation for parallel propagating linear Alfvén waves in the presence of the Cowling (Ohmic + ambipolar) diffusion is:

$$\begin{array}{c}\hfill {\omega }^2+i{k}^2{\eta }_{\mathrm{C}}\omega -k^2{v}_{\mathrm{A}}^2=0,\end{array}$$(B.1)

where ω is the frequency, k is the parallel wavenumber, and v_A is the Alfvén speed. The solution of Equation (B.1) is:

$$\begin{array}{c}\hfill \omega =\pm k{v}_{\mathrm{A}}\sqrt{1-\frac{k^2{\eta }_{\mathrm{C}}^2}{4v_{\mathrm{A}}^2}}-i\frac{k^2{\eta }_{\mathrm{C}}}{2}·\end{array}$$(B.2)

The complex part of the frequency is related to the damping owing to Cowling’s diffusion.

In this test, we considered a 1D domain in the x−direction with x ∈ [0, L] and with the same physical conditions as in the core of the prominence thread. The length of the domain is set to L = 10 km to consider a short wavelength and so to enhance the role of diffusion. The magnetic field is uniform and aligned with the x−direction as well. The fundamental standing Alfvén wave is excited by considering an initial condition for the y-component of the velocity as

$$\begin{array}{c}\hfill v_y=v_{0}sin\left(\frac{\pi x}{L}\right),\end{array}$$(B.3)

with v₀ = 0.01v_A to be in the linear regime. Regarding the boundary conditions, we set all the variables to satisfy outflow conditions, that is, zero gradient, except for the three components of velocity, which are fixed to zero. We used a numerical grid of 1000 points.

We studied the damping of the Alfvén wave by fitting the maximum value of the y−component of velocity with an exponentially damped sine function. The result is shown in Fig. B.1. We normalized all quantities using the dimensional values of L and v_A. Then, in dimensionless values k = π and η_C ≈ 1.96 × 10⁻³. Using Eq. (B.2), the theoretical frequency is ω ≈ 0.314010 − i 0.009672. In turn, the frequency obtained from the simulations is ω ≈ 0.313989 − i 0.009699. Therefore, we find that the simulations correctly recover the theoretical frequency.

Fig. B.1. Temporal evolution of the y component of velocity at the anti-node. The blue dots corresponds to the results from the 1.5D simulation while the black solid line corresponds to the fitting of the data by a damped sine function, as indicated in the plot. Code units are used.

B.2. For 2.5 dimensions

We solved the diffusion problem of the magnetic field in 2D under the physical conditions of the core of the prominence thread. The velocity components are all set to zero, so the problem is solved under static conditions. This test is inspired by that included in the PLUTO documentation to verify the implementation of the resistivity module¹.

We solved the problem in a uniform grid of 512x512 points where x, y ∈ [ − L, L], with L = 10 km. We considered a uniform background magnetic field in the z-direction with strength B₀. Outflow boundary conditions are used.

Again, we normalized all quantities using the dimensional values of L and v_A. We considered the approximate form of the resistivity tensor given in Eq. (13). The simulation is initiated at $\overline{t}=1$ with the following prescription for the components of the magnetic field:

$$\begin{array}{cc}\hfill B_x(\overline{t}=1)& =\epsilon exp(-\frac{y^2}{4{\eta }_{\mathrm{O}}}),\hfill \end{array}$$(B.4)

$$\begin{array}{cc}\hfill B_y(\overline{t}=1)& =\epsilon exp(-\frac{x^2}{4{\eta }_{\mathrm{O}}}),\hfill \end{array}$$(B.5)

$$\begin{array}{cc}\hfill B_z(\overline{t}=1)& =B_0+\epsilon exp[-\frac{x^2+y^2}{4{\eta }_C}],\hfill \end{array}$$(B.6)

with ε = 0.001B₀, so that B_x, B_y ≪ B_z to be consistent with our approximate implementation of the ambipolar diffusion. The temporal evolution of the magnetic field for $\overline{t}>1$ can analytically be obtained as

$$\begin{array}{cc}\hfill B_x(\overline{t})& =\frac{\epsilon }{\sqrt{\overline{t}}}exp(-\frac{y^2}{4{\eta }_{\mathrm{O}}\overline{t}}),\hfill \end{array}$$(B.7)

$$\begin{array}{cc}\hfill B_y(\overline{t})& =\frac{\epsilon }{\sqrt{\overline{t}}}exp(-\frac{x^2}{4{\eta }_{\mathrm{O}}\overline{t}}),\hfill \end{array}$$(B.8)

$$\begin{array}{cc}\hfill B_z(\overline{t})& =B_0+\frac{\epsilon }{\overline{t}}exp[-\frac{x^2+y^2}{4{\eta }_C\overline{t}}],\hfill \end{array}$$(B.9)

As a verification of the test, we plotted in Fig. B.2 the evolution of the z−component of the magnetic field at the center of the numerical domain, x = y = 0, obtained from the numerical simulation. According to Equation (B.9), B_z should decrease as $1/\overline{t}$ in that point. The numerical results agree perfectly with the expected dependence. We have verified that the other components of the magnetic field also follow the analytical result (not shown here), meaning that the diffusion problem of the magnetic field is correctly solved in 2D in a situation with B_x, B_y ≪ B_z.

Fig. B.2. Temporal evolution of the z−component of the magnetic field at the center of the numerical domain. The blue dots corresponds to the results from the 2.5D simulation while the red dashed line corresponds to the fitting of the data by the function: $f(\overline{t})=1.0+1.0/\overline{t}$. Code units are used.

B.3. For three dimensions

As in the Appendix B.2, we solved the diffusion problem of the magnetic field under the physical conditions of the core of the prominence thread, but now in a 3D domain. We used 128x128x128 points in a uniform grid where x, y, z ∈ [ − L, L], with L = 10 km as before. The situation is equivalent to that solved in the Appendix B.2 but adding the third dimension.

The simulation is initiated at $\overline{t}=1$ with the following prescription for the components of the magnetic field:

$$\begin{array}{cc}\hfill B_x(\overline{t}=1)& =\epsilon exp(-\frac{y^2}{4{\eta }_{\mathrm{O}}})exp(-\frac{z^2}{4{\eta }_{\mathrm{C}}}),\hfill \end{array}$$(B.10)

$$\begin{array}{cc}\hfill B_y(\overline{t}=1)& =\epsilon exp(-\frac{x^2}{4{\eta }_{\mathrm{O}}})exp(-\frac{z^2}{4{\eta }_{\mathrm{C}}}),\hfill \end{array}$$(B.11)

$$\begin{array}{cc}\hfill B_z(\overline{t}=1)& =B_0+\epsilon exp[-\frac{x^2+y^2}{4{\eta }_C}],\hfill \end{array}$$(B.12)

with ε = 0.001B₀ as before. The analytic temporal evolution of the magnetic field for $\overline{t}>1$ is:

$$\begin{array}{cc}\hfill B_x(\overline{t})& =\frac{\epsilon }{\overline{t}}exp(-\frac{y^2}{4{\eta }_{\mathrm{O}}\overline{t}})exp(-\frac{z^2}{4{\eta }_{\mathrm{C}}\overline{t}}),\hfill \end{array}$$(B.13)

$$\begin{array}{cc}\hfill B_y(\overline{t})& =\frac{\epsilon }{\overline{t}}exp(-\frac{x^2}{4{\eta }_{\mathrm{O}}\overline{t}})exp(-\frac{z^2}{4{\eta }_{\mathrm{C}}\overline{t}}),\hfill \end{array}$$(B.14)

$$\begin{array}{cc}\hfill B_z(\overline{t})& =B_0+\frac{\epsilon }{\overline{t}}exp[-\frac{x^2+y^2}{4{\eta }_C\overline{t}}],\hfill \end{array}$$(B.15)

As a verification of the 3D test, we plot in Fig. B.3 the evolution of the y−component of the magnetic field at the center of the numerical domain, x = y = z = 0. According to Equation (B.6), B_y should decrease as $1/\overline{t}$ in that point. As expected, the numerical data follows the theoretical dependence. The other components of the magnetic field are also correctly evolved (not shown here). Therefore, the diffusion problem of the magnetic field is correctly solved in 3D.

Fig. B.3. Temporal evolution of the y−component of the magnetic field at the center of the numerical domain. The blue dots corresponds to the results from the 3D simulation while the red dashed line corresponds to the fitting of the data by the function: $g(\overline{t})=1.0/\overline{t}$. Code units are used.

All Figures

	Fig. 1. Sketch of the prominence thread model. The solar photosphere is represented by the two gray planes at both ends of the magnetic flux tube.
In the text

	Fig. 2. Equilibrium density profile. Top panel: longitudinal dependence at r = 0. Bottom panel: transverse dependence at y = z = 0. We used l/R = 0.6, ρi, 0 = 100ρe, L/R = 50, and χ = 100. The density profiles are normalized with respect to the external, coronal density.
In the text

	Fig. 3. Same as Fig. 2, but for the sound and Alfvén speeds in the equilibrium. In both panels, the blue solid line corresponds to the Alfvén speed and the red dashed line corresponds to the sound speed. Both speeds are normalized with respect to the external Alfvén speed, vA, e.
In the text

	Fig. 4. Longitudinal cuts along the tube axis, r = 0, of the (a) Ohmic diffusion coefficient, (b) ambipolar diffusion coefficient, (c) Cowling diffusion coefficient, (d) temperature, and (e) ionization fraction in the equilibrium.
In the text

Fig. 5. Cross-sectional cuts of density (top left), temperature (top right), current density squared (bottom left), and vorticity squared (bottom right) at the end of the ideal MHD simulation, $\overline{t}=820$. Temperature is normalized with respect to the initial value at r = 0. Logarithmic scale is used for each variable except for the density to optimize visualization. The cross-sectional cuts are done at the tube center, z = 0, except for current density squared which is done at z = R. The complete temporal evolution is available as an online movie.

In the text

	Fig. 6. Same as Fig. 5, but the cut is done longitudinally at y = 0. Note: the horizontal and vertical axes are not to scale. The complete temporal evolution is available as an online movie.
In the text

Fig. 7. Top panel: temporal evolution of the Fourier coefficient with p = 0 (red solid line) and of the sum of the modulus of the first twenty positive Fourier coefficients (blue dashed line). For this analysis, the azimuthal component of velocity along a circle of radius 1.3R at the z = 0 plane is used. The two vertical dashed black lines show, respectively, the onset of the KHI obtained from this Fourier analysis, τKH ≈ 70, and the time at which the KHI is first seen to grow in the density, τvis = 280. Bottom panel: same as the top panel but for the radial component of velocity.

In the text

Fig. 8. Temporal evolution of the vorticity squared (top panel) and the current density squared (bottom panel) integrated in the whole computational domain. The vorticity values are normalized with respect to the initial value. The two dashed black lines show respectively the onset of the Kelvin–Helmholtz instability inferred from the analysis of azimuthal modes, τKH = 70, and the time at which the integrated vorticity squared saturates, τsat = 300.

In the text

	Fig. 9. Cross-sectional cut at z = R of the current density squared in logarithmic scale at the end of the simulations, $\overline{t}=820$. The panels are sorted as follows: ideal simulation (left), Ohmic simulation (center), and Cowling simulation (right). The complete temporal evolution is available as an online movie.
In the text

	Fig. 10. Same as bottom panel from Fig. 8 for the ideal simulation (solid red line), but now including the results from the Ohmic simulation (blue dashed line) and the Cowling simulation (black dash-dotted line).
In the text

	Fig. 11. Snapshot of the volumetric heating rate, Q, in a cross sectional cut at z = R at the end of the simulation time, $\overline{t}=820$. A logarithmic scale has been used in both panels for an optimal visualization. Left panel: Ohmic simulation. Right panel: Cowling simulation. The complete temporal evolution is available as an online movie.
In the text

	Fig. 12. Temporal evolution of the average temperature in a subdomain where x/R ∈ [−1.4,1.4], y/R ∈ [−1.4,1.4], and z/R ∈ [−4,4] for the ideal simulation (solid red line), the Ohmic simulation (blue dashed line) and the Cowling simulation (black dash-dotted line). The temperature is normalized with respect to the initial average temperature.
In the text

	Fig. 13. Normalized specific intensity at the Hα rest frequency in four different times. The normalization is done with respect to the maximum value at the beginning of the simulation. We note that the longitudinal, z, and transverse, x, directions to the thread are not to scale. A complete temporal evolution is available as an online movie.
In the text

	Fig. 14. Time-distance map of the normalized specific intensity across the thread at the rest frequency of Hα (left). The normalization is done with respect to the maximum value at the beginning of the simulation. Right panel is the same as the left, but neglecting the LOS velocities.
In the text

	Fig. A.1. Temporal evolution of the maximum value of Bx/Bz (top). The same, but for By/Bz (bottom). The results are obtained from the ideal MHD simulation.
In the text

	Fig. B.1. Temporal evolution of the y component of velocity at the anti-node. The blue dots corresponds to the results from the 1.5D simulation while the black solid line corresponds to the fitting of the data by a damped sine function, as indicated in the plot. Code units are used.
In the text

	Fig. B.2. Temporal evolution of the z−component of the magnetic field at the center of the numerical domain. The blue dots corresponds to the results from the 2.5D simulation while the red dashed line corresponds to the fitting of the data by the function: $f(\overline{t})=1.0+1.0/\overline{t}$. Code units are used.
In the text

	Fig. B.3. Temporal evolution of the y−component of the magnetic field at the center of the numerical domain. The blue dots corresponds to the results from the 3D simulation while the red dashed line corresponds to the fitting of the data by the function: $g(\overline{t})=1.0/\overline{t}$. Code units are used.
In the text