© ESO 2022

1. Introduction

Solar prominences are inhomogeneous condensations of cool plasma in the corona (e.g., Lin et al. 2005; Berger et al. 2008). The fine structure of prominences consists of very long and thin magnetic flux tubes partially filled with cold plasma (in comparison with the hot coronal plasma), which are usually referred to as threads or fibrils (Okamoto et al. 2007). The typical length of these magnetic flux tubes is of the order of 10⁵ km, with a radius between 50 km and 500 km. The length of the region of the tube filled with cold plasma inferred from Hα observations is in the range of 3000 km–28 000 km (Lin et al. 2008; Lin 2011).

Observations have revealed the presence of transverse oscillations in prominence threads, with typical periods from 1 to 20 min (Lin et al. 2007, 2009; Ning et al. 2009). These oscillations are usually interpreted as magnetohydrodynamic (MHD) kink waves (see e.g., Terradas et al. 2008b; Soler et al. 2010a; Arregui et al. 2011). Their properties have been used to perform prominence seismology (Arregui et al. 2012, 2018; Ballester et al. 2014) with the intention of inferring information about the physical parameters of the plasma and/or the magnetic field (Montes-Solís & Arregui 2019).

For instance, the comparison of the periods of the fundamental mode, P₀, and the corresponding first overtone, P₁, can be used to obtain information about the longitudinal distribution of density. The ratio of the damping time to the period has commonly been used to infer the radial distribution and determine which physical mechanism is responsible for that attenuation. In a longitudinally homogeneous tube, the period ratio is P₀/P₁ = 2 (Edwin & Roberts 1983), but it differs from 2 if a longitudinal inhomogeneity is considered. Andries et al. (2005a,b) showed that for the case of coronal loops, which are denser at the footpoints than at the apex, the period ratio is smaller than 2. Later, Díaz et al. (2010) showed that for a prominence thread, which is denser at the centre of the tube, the period ratio is always greater than 2, a result that was confirmed by Arregui et al. (2011) and Soler & Goossens (2011), for example. The relation P₀/P₁ has been used to show that the ratio between the average density of the tube and the density at its centre is very small, which implies that there are very strong gradients of density along the threads (Soler et al. 2015). Regarding the damping-time-to-period ratio, Soler et al. (2009, 2014) analysed several possible damping mechanisms and determined that the process that better describes the data from observations of prominence threads is resonant absorption (Ionson 1978; Ruderman & Roberts 2002; Arregui et al. 2008; Soler et al. 2010a). This phenomenon occurs when there is a smooth variation between the internal and external densities of the flux tube and depends on the thickness of this radial transition layer and on the density contrast between the internal and external plasmas.

The previously mentioned works usually focused on the density inhomogeneity of the threads in one direction only, namely the longitudinal or the radial one. In that way, analytical expressions could be provided for the periods and damping times of oscillations in tubes with longitudinal uniform density or represented by piecewise constant models. On the contrary, the consideration of the inhomogeneity along those two directions at the same time requires a numerical treatment (Arregui et al. 2005, 2011).

In the present investigation, we used a model of a prominence thread in which the density follows a Lorentzian profile along the longitudinal direction. The main parameter that describes this profile is the longitudinal density ratio, χ, which is defined as the ratio between the densities at the centre and at the ends of the tube. For prominence threads, this parameter takes values larger than 1. In the radial direction there is a smooth transition layer that connects the internal plasma with the external plasma. The variation of density in this radial transition layer is given by a sinusoidal profile. Then, we studied the properties of transverse kink oscillations by applying two different methods. Firstly, we used the numerical code MoLMHD (Terradas et al. 2016, 2018), which solves the non-linear ideal MHD equations, to run 3D simulations and compute the full temporal evolution of the flux tube. Then, we compare the obtained results with the solutions from a 2D generalised eigenvalue problem (Soler et al. 2015). Numerical simulations are very computationally costly, but they have the advantage that they capture the non-linear dynamics of this kind of oscillation. In contrast, the eigenvalue problem is restricted to the linear regime, but it is much faster and it is more appropriate for performing a parametric analysis of the periods and damping times. Therefore, both methods complement each other.

We studied how the periods, damping times, and amplitudes of the fundamental and first harmonic modes vary with the longitudinal density ratio and with the average density. We computed the period ratio P₀/P₁ and the damping to period ratio, and we compared the obtained results with those corresponding to tubes with longitudinally uniform profiles of density.

After this analysis, we proceeded in our numerical research and used the results from the 3D simulations to perform forward modelling and provide synthetic profiles in the Hα line. We used the approximate method described in Heinzel et al. (2015). We computed the Hα intensity along two different lines of sight (parallel and perpendicular to the direction of oscillation) and provide snapshots that show the longitudinal and radial distribution of the density of prominence threads. In this paper, we discuss the observational signatures of the transverse oscillations and how the strong variation of density along the thread may have implications for prominence seismology.

The synthetic Hα profiles do not only provide information about the amplitudes, periods and damping times of the kink oscillations. They also include evidences of non-linear processes that are taking place in the threads. For instance, large velocity gradients at the boundaries between the internal and external layers of the flux tube trigger shear instabilities, such as the Kelvin-Helmholtz instability (KHI, Chandrasekhar 1961; Browning & Priest 1984; Terradas et al. 2008a). As a consequence, large deformations appear at those boundaries as the instability develops (see e.g., Antolin et al. 2014; Magyar et al. 2015; Terradas et al. 2018). We discuss how the variations in intensity of the Hα profiles are related to these non-linear effects.

The present work is organised as follows. In Sect. 2, we describe the model we used to represent an inhomogeneous flux tube, the equations that govern its dynamics, and the methods applied to solve those equations. In Sect. 3, we analyse the characteristics of the transverse oscillations, comparing the cases of uniform and non-uniform tubes. Then, in Sect. 4, we describe the forward modelling performed to compute synthetic Hα profiles and show the observational signatures of several linear and non-linear features of the oscillations. Finally, our conclusions and a brief discussion on possible improvements are presented in Sect. 5.

2. Model and numerical setup

2.1. Flux tube model

We represent prominence threads as straight flux tubes with inhomogeneous density embedded in a hotter and lighter background (which represents the solar corona). In the longitudinal direction, the density follows a Lorentzian profile:

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

where ρ_i, 0 is the density at z = 0, L_z is the length of the tube, and χ = ρ_i, 0/ρ_i(z = L_z/2) is the longitudinal density ratio. The average density of this Lorentzian profile is given by

$$\begin{array}{c}\hfill ⟨{\rho }_{\mathrm{i}}⟩=\frac{1}{L_z}{\displaystyle {\int }_{-L_z/2}^{L_z/2}{\rho }_{\mathrm{i}}(z)\mathrm{d}z={\rho }_{\mathrm{i},0}\frac{arctan\sqrt{\chi -1}}{\sqrt{\chi -1}}.}\end{array}$$(2)

According to this formula, the value of ⟨ρ_i⟩ decreases as χ increases.

In the transverse direction to the tube, the density is radially uniform inside the tube with a smooth transition between the internal and external values of density, ρ_in and ρ_ex, respectively. In this transition layer, the density varies as

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

with $r=\sqrt{x^2+{\mathit{y}}^2}$, R being the radius of the tube, and l being the width of the transition layer. In this work, the relations between the internal and external densities and between the parameters l and R are set to ρ_i, 0 = 100ρ_ex (where ρ_ex = ρ_ref and ρ_ref is the reference value of the coronal density) and l/R = 0.3, respectively.

The flux tube model is displayed in Figs. 1 and 2. A 3D sketch of an inhomogeneous thread with χ = 10 is shown in the former: the red, green, and blue surfaces correspond to different density contours, and the vertical grey lines show the magnetic field lines. Figure 2 displays the vertical distribution of density for different values of the density ratio (top panel) and its radial profile at the tube’s mid-height (bottom panel).

Fig. 1. Sketch of equilibrium state for a thread with a density ratio of χ = 10. Density contours for ρ = 25ρex (red), ρ = 50ρex (green), ρ = 75ρex (dark blue), and ρ = 100ρex (light blue). Magnetic field lines, with line-tying conditions applied at the top and bottom planes, are shown in grey.

Fig. 2. Density profiles of equilibrium state. Top panel: vertical profiles for several values of the density ratio: χ = 5 (black), χ = 10 (green), χ = 20 (light blue), χ = 50 (dark blue), and χ = 100 (red). Bottom panel: radial profile at z = 0.

To perform the simulations, we considered a numerical domain of −5 < x/L₀ < 5, −5 < y/L₀ < 5, and −L_z/2 < z/L₀ < L_z/2, where L₀ = 1 Mm is the reference length of the system and L_z = 50R, with R = L₀. We note that the considered ratio L_z/R produces a tube that is thicker than what is expected for a prominence thread. For instance, the observations performed by Okamoto et al. (2007), Lin et al. (2008), and Lin (2011) provided much larger values of the ratio L_z/R, which range from 400 to 2000. The reason why we chose this small ratio is related to the computing time of the simulations: the oscillation period of a flux tube is proportional to its length, as shown by Edwin & Roberts (1983), so the study of transverse oscillations of tubes with much larger values of L_z/R requires much longer computing times. The investigation we present in this paper required a considerable number of simulations, so it would have been impractical to use tubes with much higher values of L_z/R. Nevertheless, we expect that this choice does not affect the general validity of our results, as they mainly depend on the longitudinal variation of the tube properties and not on its actual length or radius. In addition, we considered the presence of a uniform longitudinal magnetic field of B₀ = 10 G, which permeates the whole domain, and, in combination with the reference value of the coronal density ρ_ref = 3.656 × 10⁻¹³ kg m⁻³, it gives a reference Alfvén speed of $c_{\mathrm{A},\mathrm{ref}}=B_0/\sqrt{\mu {\rho }_{\mathrm{ref}}}\approx 1494\mathrm{km}{\mathrm{s}}^{-1}$ and a reference time t₀ = L₀/c_A,ref ≈ 0.67 s.

The pressure is uniform in the whole domain, which implies that the system is initially in mechanical equilibrium. We used a value of P = 0.005 Pa, which corresponds to a plasma beta of $\beta =2{\mu }_{0}P/B_0^2\approx 0.013$. In this low-β limit, pressure has a negligible effect on the linear regime of transverse oscillations. However, it becomes relevant in the non-linear regime. Hollweg (1971), Rankin et al. (1994, 1995), Tikhonchuk et al. (1995), and Terradas & Ofman (2004) demonstrated that the ponderomotive force causes longitudinal flows that tend to accumulate the mass at the centre of the tube. In a pressureless plasma, this accumulation would go on without bound. On the contrary, the effect of pressure is to prevent this unlimited growth and turn it into an oscillatory behaviour.

We applied line-tying boundary conditions at the ends of the flux tube to represent the field lines anchoring in the photosphere, meaning that the three components of the velocity are set to zero and the normal component of the magnetic field is fixed, while the normal derivatives of the rest of variables are equal to zero. For the lateral planes, closed boundary conditions are chosen. We note that these conditions may produce unwanted reflections of perturbations, but we set the lateral boundaries far enough from the numerical domain of interest so they have a negligible influence on the evolution of the flux tube.

2.2. Methods and equations

On the one hand, we performed 3D simulations using the numerical code MoLMHD (Terradas et al. 2016, 2018), which solves the non-linear ideal MHD equations:

$$\begin{array}{c}\hfill \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·\left(\rho \mathit{V}\right)=0,\end{array}$$(4)

$$\begin{array}{c}\hfill \frac{\partial \left(\rho \mathit{V}\right)}{\partial t}+\mathrm{\nabla }·(\rho \mathit{VV}+P\mathbb{I}-\frac{\mathit{BB}}{\mu }+\frac{{\mathit{B}}^2}{2\mu }\mathbb{I})=0,\end{array}$$(5)

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}=\mathrm{\nabla }\times (\mathit{V}\times \mathit{B}),\end{array}$$(6)

$$\begin{array}{c}\hfill \frac{\partial P}{\partial t}+\mathrm{\nabla }·\left(\gamma P\mathit{V}\right)=(\gamma -1)\mathit{V}·\mathrm{\nabla }P,\end{array}$$(7)

where ρ, V, B, and P represent the density, velocity, magnetic field, and pressure, respectively; 𝕀 is the unit tensor; μ is the magnetic permeability; and γ is the adiabatic index. These equations are discretised using a combination of a sixth-order central finite-differences scheme and a weighted essentially non-oscillatory (WENO) scheme (see details in Terradas et al. 2015) and are advanced in time through a third-order total variation diminishing (TVD) Runge-Kutta method (Harten 1983).

On the other hand, we solved an updated version of the 2D generalised eigenvalue problem described in Soler et al. (2015). That work did not take into account the presence of the transition layer between the internal and external densities but considered the existence of a density jump. Therefore, it could not predict the damping of the oscillations associated with the phenomenon of resonant absorption. For the present work, the effect of transverse inhomogeneity has been included in the solution of the eigenvalue problem. For this method, we linearised the basic MHD Eqs. (4)–(7). We only considered transverse oscillations, which are accurately described in the β = 0 approximation. Therefore, we drop the effect of gas pressure for simplicity. Unlike in Soler et al. (2015), the considered model is transversely non-uniform, and the ideal global mode perturbations are singular (resonant) at the specific radial position where the global mode frequency matches the local Alfvén frequency. In order to avoid the ideal singularity, which is not tractable numerically, here we need to add a small diffusive term in the induction equation. This term does not alter the overall behaviour of the perturbations, nor does it modify the eigenfrequencies, as long as it is sufficiently small (see e.g., Arregui et al. 2011). So, the relevant linearised equations for the computation of eigenvalues are

$$\begin{array}{c}\hfill {\rho }_0\frac{\partial \mathit{v}}{\partial t}=\frac{1}{\mu }(\mathrm{\nabla }\times \mathit{b})\times {\mathit{B}}_{\mathbf{0}},\end{array}$$(8)

$$\begin{array}{c}\hfill \frac{\partial \mathit{b}}{\partial t}=\mathrm{\nabla }\times (\mathit{v}\times {\mathit{B}}_{\mathbf{0}})-\eta \mathrm{\nabla }\times \mathrm{\nabla }\times \mathit{b},\end{array}$$(9)

where v = (v_r,v_ϕ,0) and b = (b_r,b_ϕ,b_z) are the velocity and magnetic field perturbations, respectively, expressed in a cylindrical coordinate system with the z-axis aligned with the flux tube axis, ρ₀ and B₀ = (0, 0, B₀) are the background density and magnetic field, respectively, and η is the coefficient of magnetic diffusion (assumed uniform). We put the perturbations proportional to exp(imϕ), where m is the azimuthal wave number, and set m = 1 corresponding to kink modes. In addition, we express the temporal dependence of the perturbations as exp(−iωt), where ω is the complex eigenfrequency. The period of the oscillation is P = 2π/Re(ω) and the exponential damping time is τ_D = −1/Im(ω). Then, Eqs. (8) and (9) are expanded as

$$\begin{array}{c}\hfill \omega v_r=i\frac{B_0}{\mu {\rho }_0}(\frac{\partial b_r}{\partial z}-\frac{\partial b_z}{\partial r}),\end{array}$$(10)

$$\begin{array}{c}\hfill \omega v_{\varphi }=\frac{B_0}{\mu {\rho }_0}(\frac{b_z}{r}+i\frac{\partial b_{\varphi }}{\partial z}),\end{array}$$(11)

$$\begin{array}{cc}& \omega b_r=i{B}_0\frac{\partial v_r}{\partial z}-i\eta (\frac{b_r}{r^2}+i\frac{1}{r}\frac{\partial b_{\varphi }}{\partial r}+i\frac{b_{\varphi }}{r^2}\hfill \\ \hfill & -\frac{{\partial }^2{b}_r}{\partial z^2}+\frac{{\partial }^2{b}_z}{\partial r\partial z}),\hfill \end{array}$$(12)

$$\begin{array}{cc}& \omega b_{\varphi }=i{B}_0\frac{\partial v_{\varphi }}{\partial z}+i\eta (\frac{{\partial }^2{b}_{\varphi }}{\partial r^2}+\frac{1}{r}\frac{\partial b_{\varphi }}{\partial r}-\frac{b_{\varphi }}{r^2}\hfill \\ \hfill & -i\frac{1}{r}\frac{\partial b_r}{\partial r}+i\frac{b_r}{r^2}-i\frac{1}{r}\frac{\partial b_z}{\partial z}+\frac{{\partial }^2{b}_{\varphi }}{\partial z^2}),\hfill \end{array}$$(13)

$$\begin{array}{cc}& \omega b_z=-i{B}_0(\frac{\partial v_r}{\partial r}+\frac{v_r}{r}+i\frac{v_{\varphi }}{r})+i\eta (\frac{{\partial }^2{b}_z}{\partial r^2}+\frac{1}{r}\frac{\partial b_z}{\partial r}\hfill \\ \hfill & -\frac{b_z}{r^2}-\frac{{\partial }^2\partial b_r}{\partial r\partial z}-\frac{1}{r}\frac{\partial b_r}{\partial z}-i\frac{1}{r}\frac{\partial b_{\varphi }}{\partial z}).\hfill \end{array}$$(14)

Following Soler et al. (2015), trapped kink oscillations must satisfy the following boundary conditions:

$$\begin{array}{c}\hfill \frac{\partial v_r}{\partial r}=\frac{\partial v_{\varphi }}{\partial r}=\frac{\partial b_r}{\partial r}=\frac{\partial b_{\varphi }}{\partial r}=0,b_z=0,\mathrm{at}r=0,\end{array}$$(15)

$$\begin{array}{c}\hfill v_r=v_{\varphi }=b_r=b_{\varphi }=b_z=0,\mathrm{at}r=r_{\max },\end{array}$$(16)

$$\begin{array}{c}\hfill v_r=v_{\varphi }=b_z=0,\frac{\partial b_r}{\partial z}=\frac{\partial b_{\varphi }}{\partial z}=0,\mathrm{at}z=\pm L_z/2,\end{array}$$(17)

where r = 0 is the centre of the tube and r = r_max is a radial position sufficiently far away from the tube. Ideally, r_max → ∞, while here we considered r_max = 20R for practical purposes. Equations (10)–(14) together with the above boundary conditions define a 2D generalised eigenvalue problem, where ω is the eigenvalue and the perturbations form the vector of eigenfunctions. The eigenvalue problem is solved with a more efficient numerical method than that used in Soler et al. (2015). The scheme is implemented in a Wolfram Mathematica routine, where the spatial integration of the equations in the r and z directions is done with finite elements in a 2D structured mesh with a resolution of 400 × 40 cells. The mesh resolution is non-uniform in the radial direction, as a much finer resolution is needed in the non-uniform transitional layer to correctly describe the spatial scales associated with the resonant behaviour of the perturbations. The routine returns the closest eigenvalue to an initially provided guess and can be run iteratively to compute the solution as a function of a specific parameter of the model.

2.3. Initial perturbation

Transverse kink oscillations in flux tubes have the azimuthal wave number m = 1 (Edwin & Roberts 1983; Nakariakov & Verwichte 2005). To trigger this kind of oscillations in our simulations, we applied the following perturbation to the x-component of velocity:

$$\begin{array}{c}\hfill V_x(t=0)=V_{0}f(r)cos(n{k}_{z}z+\frac{(n-1)}{2}\pi ),\end{array}$$(18)

where

$$\begin{array}{c}\hfill f(r)=\frac{\rho (r,z=0)-{\rho }_{\mathrm{ex}}}{{\rho }_{\mathrm{i},0}-{\rho }_{\mathrm{ex}}},\end{array}$$(19)

k_z = π/L_z and n is the longitudinal oscillation mode. The amplitude of the perturbation is expressed in terms of the internal Alfvén speed, V₀ = αc_A, in, where $c_{\mathrm{A},\mathrm{in}}=B_0/\sqrt{\mu {\rho }_{i,0}}$ and α is a dimensionless factor that controls the strength of the perturbation. The factor f(r) is used to ensure that the perturbation is applied to the tube and the transition layer but not to the external background plasma.

If n = 1, the perturbation corresponds to the fundamental longitudinal mode of oscillation; n = 2 corresponds to the first longitudinal overtone. In a tube with a uniform density, the periods of these modes in the thin tube approximation are given by (Edwin & Roberts 1983)

$$\begin{array}{c}\hfill P_0=\frac{2L_z}{c_{\mathrm{A},\mathrm{in}}}\sqrt{\frac{1+{\rho }_{\mathrm{ex}}/{\rho }_{\mathrm{i},0}}{2}}\mathrm{and}{P}_1=\frac{P_0}{2},\end{array}$$(20)

respectively. In the following sections, we check how the inhomogeneity in density modifies the oscillation periods of these two modes.

We note that the longitudinal dependence given in Eq. (18) corresponds to an eigenfunction of kink oscillations in longitudinally uniform tubes but not for the case with a longitudinal inhomogeneity (Soler et al. 2010a, 2015; Arregui et al. 2011). Therefore, this perturbation will trigger additional longitudinal harmonics. However, we expect the contribution of these modes to have a negligible effect on the evolution of the transverse oscillations we are interested in analysing.

3. Analysis of numerical results

In this section, we present the results of our numerical study of oscillations of inhomogeneous threads and provide comparisons with the case of longitudinally homogeneous threads. All the simulations analysed here used a numerical domain of 200 × 200 × 100 points, which corresponds to a spatial resolution of 50 km in the x and y directions, and of 500 km in the vertical direction.

3.1. Example simulations

We start by describing a simulation of the fundamental longitudinal mode of oscillation in an inhomogeneous thread with a longitudinal density ratio of χ = 10. This transverse oscillation is triggered by applying the initial perturbation given by Eq. (18) with n = 1. For this example, the amplitude of the perturbation is V₀ = 0.05c_A, in.

The perturbation causes a lateral displacement of the tube that follows a sinusoidal motion. The displacement is maximum at the central position of the tube and zero at its ends. This behaviour is illustrated by the density snapshots displayed in Fig. 3. In addition to the oscillation around the equilibrium position, the horizontal cuts of density at the height z = 0 (top panels of Fig. 3) show that the shape of the originally cylindrical thread is deformed as time advances. This is a consequence of the triggering of the KHI (Terradas et al. 2008a; Antolin et al. 2014; Magyar et al. 2015) and the non-linear generation of fluting modes (Ruderman & Goossens 2014; Magyar & Van Doorsselaere 2016; Terradas et al. 2018), which are modes with the azimuthal wave number m ≥ 2 (Nakariakov & Verwichte 2005) localised around the boundary of the tube (Soler 2017).

Fig. 3. Snapshots of a simulation with χ = 10, V0 = 0.05cA, in, n = 1, and a numerical domain of 200 × 200 × 100 points. Top: density colour maps at z = 0. Bottom: density colour maps at y = 0. Dashed lines represent the density contour ρ = 50ρref at t = 0.

Now, we turn our attention to the case of the first longitudinal overtone, with n = 2. For this simulation, we used the same value of V₀ as in the previous case. Figure 4 shows that the first overtone has a node at z = 0. There is no lateral displacement at the position where the tube is denser. Conversely, it has two anti-nodes, one at each half of the tube, which oscillate laterally in opposing directions. We see that the amplitude of the displacement for this mode is smaller than that for the fundamental mode although the same value of V₀ was used, and, consequently, the boundary of the tube is less distorted.

Fig. 4. Snapshots of a simulation with χ = 10 and n = 2. Top: density colour maps at z = −Lz/4 = −12.5L0. Bottom: density colour maps at y = 0. Dashed lines represent the density contours ρ ≈ 15ρref (top panels) and ρ = 50ρref (bottom) at t = 0.

In order to obtain a more precise description of the two simulated oscillations, we computed the displacement of the centre of mass of the density cuts located at z = 0 and z = −L_z/4 = −12.5L₀, respectively. The results are presented in Fig. 5. The top panel confirms that for the mode n = 1, the different longitudinal positions of the tube oscillate in phase, and it shows that the amplitude of the oscillation decreases with time. The main reason for this attenuation is the linear process of resonant absorption, which causes a transfer of energy from the global kink oscillation to localised Alfvén modes in the transverse inhomogeneous layer (Ionson 1978; Ruderman & Roberts 2002; Arregui et al. 2008; Soler et al. 2010a). The non-linear excitation of fluting modes (see e.g., Ruderman & Goossens 2014; Soler 2017) and the triggering of the KHI (Terradas et al. 2018; Van Doorsselaere et al. 2021) also contribute to the damping of the oscillation.

Fig. 5. Displacement of tube axis at different longitudinal positions from simulations of modes n = 1 (top) and n = 2 (bottom). The solid and dotted lines represent oscillations with periods given by Eq. (20) for tubes of uniform density equal to ρi, 0 and ⟨ρi⟩, respectively. The dashed lines represent the exponential decay of the displacement of the fundamental mode.

The bottom panel of Fig. 5 shows that for the first overtone there is no lateral displacement at z = 0. The displacement at z = −L_z/4 follows a similar behaviour than that of the fundamental mode but with smaller amplitude and shorter period.

In Fig. 5, we include black solid lines that represent oscillations of periods P₀ and P₁, respectively, which are given by Eq. (20). These periods correspond to kink oscillations in longitudinally uniform thin tubes with density equal to ρ_i, 0. The dotted black lines represent oscillations of tubes with longitudinally uniform density given by the average density of the Lorentzian profile, that is, ρ_i(z) = ⟨ρ_i⟩. We see that the thread with χ = 10 oscillates faster than the tube with ρ_i(z) = ρ_i, 0, but slower than the tube with ρ_i(z) = ⟨ρ_i⟩, in agreement with the findings of Soler et al. (2015). The oscillation periods of the tube with the Lorentzian profile are P(n = 1)≈568t₀ and P(n = 2)≈226t₀. Therefore, the period ratio follows the relation P(n = 1)/P(n = 2) > 2, in agreement with the findings of Díaz et al. (2010), Arregui et al. (2011), and Soler & Goossens (2011) for prominence threads.

The dashed lines in the top panel of Fig. 5 show that the attenuation is described well by an exponential decay during the whole course of the simulation. This behaviour seems to differ from the results for longitudinally uniform threads obtained by Pascoe et al. (2012), Ruderman & Terradas (2013), Magyar & Van Doorsselaere (2016), and Terradas et al. (2018), who found that during the first periods of the oscillation the damping is better fitted by a Gaussian profile. Ruderman & Terradas (2013) and Magyar & Van Doorsselaere (2016) found that the damping changes from a Gaussian to an exponential profile at a switch time given by

$$\begin{array}{c}\hfill t_{\mathrm{s}}=\frac{{\rho }_{\mathrm{i},0}/{\rho }_{\mathrm{ex}}+1}{{\rho }_{\mathrm{i},0}/{\rho }_{\mathrm{ex}}-1}P.\end{array}$$(21)

The switch time decreases with the density contrast between the internal and external plasmas. For the simulations we analyse here (with ρ_i, 0 = 100ρ_ex), the switch time is approximately equal to one period. Therefore, the behaviour represented in the top panel of Fig. 5 may agree with the findings for longitudinally uniform tubes, and the apparent discrepancy may be a consequence of the chosen value for the ratio ρ_i, 0/ρ_ex instead of a consequence of the longitudinal variation of density given by the parameter χ.

3.2. Parametric study

In the following sections, we vary the value of the parameter χ (which implies a variation of the average density of the threads) and study the general behaviour of the periods, damping times, and amplitudes of the transverse oscillations.

3.2.1. Fundamental longitudinal mode

The top left panel of Fig. 6 presents the results from a series of simulations in which the initial perturbation corresponds to the fundamental mode of oscillation, n = 1, and its amplitude is given by V₀ = 0.05c_A, in. It shows the transverse displacement of the centre of mass of the tube as a function of time for the density profiles plotted in Fig. 2; black crosses, green asterisks, light blue diamonds, dark blue triangles, and red squares correspond to the longitudinal density ratios χ = 5, χ = 10, χ = 20, χ = 50, and χ = 100, respectively. We see in this panel that the amplitude, period, and the damping time of oscillation depend on the chosen longitudinal density ratio. To extract these oscillation parameters from the data of the simulations, we fitted each line to an exponentially decaying sinusoidal function:

$$\begin{array}{c}\hfill \frac{\xi }{R}=A_{0}sin\left(\frac{2\pi }{P}t\right)exp(-\frac{t}{{\tau }_{\mathrm{D}}}),\end{array}$$(22)

Fig. 6. Results of the parametric study of transverse oscillations. Top left: displacement of centre of mass as a function of time from simulations with n = 1 and with χ = 5 (black crosses), χ = 10 (green asterisks), χ = 20 (light blue diamonds), χ = 50 (dark blue triangles), and χ = 100 (red squares). Top right: oscillation period as a function of the density contrast, χ. Bottom left: damping time as a function of χ. Bottom right: ratio of damping time to period as a function of χ. In the top left panel, solid lines represent fits of the numerical results. In the remaining panels, the black solid and dashed lines correspond to the results obtained by solving the 2D eigenvalue problem for n = 1 and n = 2, respectively. Symbols represent the results from the simulations. For n = 1, Pref = P0; for n = 2, Pref = P1.

where A₀ ≡ (ξ/R)₀ is the dimensionless amplitude of the oscillation, P is the period, and τ_D is the exponential damping time. Due to non-linear effects related to the initial large amplitudes of the perturbations, some of the oscillations would be better fitted by a Gaussian damping profile, as shown by Magyar & Van Doorsselaere (2016). Nevertheless, we applied the exponential fitting function to every oscillation to be able to compare the results from the simulations with the solutions from the eigenvalue problem and with results from previous studies.

The red asterisks in the top right, bottom left, and bottom right panels of Fig. 6 display the period, damping time, and damping ratio, respectively, as functions of the longitudinal density ratio, χ. We see in the top right panel that greater density ratios correspond to shorter periods, with the case of χ = 100 having a period that is ∼0.55 times the period of the uniform tube (χ = 1). This is the same behaviour as the one presented in Fig. 3a of Soler et al. (2015). The bottom left panel shows that the damping times also decrease as the density ratio is increased. Finally, the bottom right panel shows that the dependence of the ratio between the damping time and the period on the longitudinal density ratio is not so strong as for the previous parameters, but a small increase is found as χ is increased.

The above described behaviour of the damping-time-to-period ratio is different from that of threads with no density variation along the axis. Ruderman & Roberts (2002), Arregui et al. (2008), and Soler et al. (2009, 2010a) found that in a longitudinally homogeneous thread with a radial transition layer described by Eq. (3), the ratio of damping time to period of the kink oscillations in the thin tube and thin boundary regime is

$$\begin{array}{c}\hfill \frac{{\tau }_{\mathrm{D}}}{P}\approx \frac{2}{\pi }\frac{1}{l/R}\left(\frac{{\rho }_{\mathrm{i},0}+{\rho }_{\mathrm{ex}}}{{\rho }_{\mathrm{i},0}-{\rho }_{\mathrm{ex}}}\right).\end{array}$$(23)

A similar result was obtained by Andries et al. (2005b) and Arregui et al. (2005) in their studies of coronal loop oscillations, although they considered the effect of longitudinal density stratification. The reason for this different behaviour is that in the present work the density contrast between the internal and external layers varies with the parameter χ, while in the above-mentioned studies the density contrast was a constant (Dymova & Ruderman 2006).

The model we present in Sect. 2.1 is more closely related to that used by Arregui et al. (2011), who represented a prominence thread as a tube divided in three longitudinally uniform regions with different values of density connected by small transition layers. Arregui et al. (2011) found that the damping ratio slightly increases for very short threads. In our model, short threads are represented by large values of χ.

Using Eq. (23) as a reference, the slight increase of τ_D/P with χ can be roughly understood as follows. As shown by Eq. (2), a larger χ corresponds to a smaller average density. If we use the variable ⟨ρ_i⟩ instead of ρ_i, 0 in Eq. (23), we find that τ_D/P increases as the average density decreases.

The solid lines in the top right and the bottom panels of Fig. 6 show the results from the 2D eigenvalue problem. There is a remarkable agreement with the results from simulations for the case of the oscillation periods. We see some differences in the damping times, with the simulations yielding smaller values, which means that the oscillations are attenuated slightly faster. As mentioned in Sect. 3.1, this could be caused by the non-linear generation of fluting modes (Ruderman et al. 2010) and the triggering of the KHI Terradas et al. (2018), Van Doorsselaere et al. (2021), which grows by extracting energy from the flow. The numerical diffusivity used in the code to prevent the appearance of numerical noise also adds another damping mechanism. None of these processes (either physical or numerical) are taken into account by the 2D eigenvalue problem.

3.2.2. First longitudinal harmonic

In this section, we use the same amplitude of the velocity perturbation as in the previous one, but here the perturbation corresponds to the oscillation mode n = 2, that is, the first longitudinal harmonic. This oscillation mode presents a node at the middle position of the tube, so its centre of mass remains static and cannot be used to analyse the displacement of the tube. For this mode, the maximum displacements occur at z = −L_z/4 and z = L_z/4. We chose the former position for the present analysis. We computed the centre of mass contained in the horizontal plane at that position and applied Eq. (22) to fit its displacement. The results from this fitting procedure are represented by blue diamonds on the top right and bottom panels of Fig. 6.

As for the fundamental mode, both the oscillation periods and the damping times reduce as the longitudinal density ratio is increased. However, in this case the dependence on that parameter is even more pronounced than for the case of the fundamental mode. This behaviour qualitatively agrees with the results of Andries et al. (2005a), who found that the first longitudinal overtone is more affected by the longitudinal variations of density than the fundamental mode. Then, Fig. 6 shows that the period when χ = 100 is ∼0.3 times the oscillation period of the tube with uniform density, and the damping time is reduced to less than half. This means that the inhomogeneity in density has a greater effect on the periods than on the damping times. In addition, the ratio of damping time to period has a larger increase with χ than for the case with n = 1. This is related to the large decrease of density contrast between the internal and external plasmas at z = −L_z/4 as the parameter χ increases.

3.2.3. P₀/P₁ ratio

We compared the periods of the two oscillation modes previously analysed. The results are shown in Fig. 7. For the case of a tube with uniform density, this ratio is P₀/P₁ = 2 in the thin tube approximation (Edwin & Roberts 1983). However, the value of the ratio varies when longitudinal inhomogeneities in density are taken into account. For the case of a prominence thread, the ratio is always larger than 2 (Díaz et al. 2010). The left panel of Fig. 7 shows that the period ratio grows as χ increases.

Fig. 7. Period ratio as a function of density ratio (left) and as a function of the average density of the tube (right). Dashed lines correspond to the case of tubes with uniform density. Dotted lines represent the approximation given by Eq. (24).

Up to this point, we focused on the dependence of the oscillation parameters on the ratio of the central density to the footpoint density. However, we note that each one of the profiles shown in Fig. 2 represents a thread with a different total mass, so the comparison between each other and with the uniform density case is not straightforward. A more appropriate procedure would be to consider the dependence on the total mass of the tube or, equivalently, on its average density (see Andries et al. 2005a).

The dependence of the period ratio on the average density is shown on the right panel of Fig. 7. For a longitudinally uniform thread, the period ratio is independent of the average density. For a tube with a Lorentzian profile, it increases as the average density decreases. As shown by Soler et al. (2015), this behaviour is not exclusive to the Lorentzian profile, but it is shared by other non-uniform density distributions.

The analytical approximation displayed on Fig. 7 is given by

$$\begin{array}{c}\hfill \frac{P(n=1)}{P(n=2)}\approx 1+{\left(\frac{4}{{\pi }^2}\chi \right)}^{1/4},\end{array}$$(24)

which was derived by Soler et al. (2015) after obtaining an empirical fit as a function of the average density,

$$\begin{array}{c}\hfill \frac{P(n=1)}{P(n=2)}\approx 1+{\left(\frac{⟨{\rho }_{\mathrm{i}}⟩}{{\rho }_{\mathrm{i},0}}\right)}^{-1/2},\end{array}$$(25)

using the relation given by Eq. (2) and applying the limit χ ≫ 1. Analytical proof of the validity of this formula was provided by Ruderman et al. (2016).

We see that the numerical results have a very good agreement with the analytical approximation for large values of the density ratio (or small average densities), while small differences appear in the opposite range. Prominence threads are expected to have high density ratios, so the use of this approximate formula to perform prominence seismology is justified (see e.g., Arregui & Soler 2015).

3.2.4. τ_D, 0/τ_D, 1 ratio

According to Eq. (23), the damping-time-to-period ratio in a longitudinally homogeneous tube is independent of the oscillation mode n as long as the longitudinal wavelength remains much longer than the radius of the tube, that is, in the thin tube limit. Therefore, the ratio between the damping times of the fundamental mode and the first harmonic is the same as the ratio between the periods, that is, equal to 2. However, it has been shown that when the tube has a Lorentzian profile in density, the period ratio approximately follows the relations given by Eqs. (24) and (25). We expect that these formulas are not applicable to the ratio of damping times, since the attenuation process is mainly affected by the transverse variation of density (see e.g., Ruderman & Roberts 2002; Arregui et al. 2008; Soler et al. 2009), while the oscillation period is weakly affected by this variation but strongly depends on the longitudinal one (Andries et al. 2005b; Díaz et al. 2010; Soler et al. 2015).

Figure 8 shows the damping times ratio as a function of χ (left) and on the average density (right). In addition to the results from the numerical simulations (red asterisks), this plot includes the solutions of the 2D eigenvalue problem (black solid line). We see that the damping ratio rises as the longitudinal density ratio is increased (or the average density is reduced) but not as fast as the period ratio (red dotted lines).

Fig. 8. Damping time ratio as a function of density ratio (left) and as a function of the average density of the tube (right). Red asterisks represent the results from simulations. Solid black lines correspond to the results from the eigenvalue problem. Horizontal dashed lines represent the case of tubes with uniform density. The red dotted lines represent the approximations given by Eqs. (24) and (25). The blue dotted-dashed lines in the left and right panels represent the fits given by Eqs. (26) and (27), respectively.

With the goal of finding expressions similar to Eqs. (24) and (25) for the damping ratio, we fitted the solutions of the eigenvalue problem displayed in Fig. 8. We used the fitting function f(X) = 1 + a₁X^a₂, where f(X) is the ratio of the damping times, X is the independent variable (χ or the ratio ⟨ρ_i⟩/ρ_i, 0), and a₁ and a₂ are the coefficients of the fit. We obtained the following approximations:

$$\begin{array}{c}\hfill \frac{{\tau }_{\mathrm{D}}(n=1)}{{\tau }_{\mathrm{D}}(n=2)}\approx 1+{\chi }^{0.15}\end{array}$$(26)

and

$$\begin{array}{c}\hfill \frac{{\tau }_{\mathrm{D}}(n=1)}{{\tau }_{\mathrm{D}}(n=2)}\approx 1+1.05{\left(\frac{⟨{\rho }_{\mathrm{i}}⟩}{{\rho }_{\mathrm{i},0}}\right)}^{-0.34},\end{array}$$(27)

which are shown in Fig. 8 as light blue dotted-dashed lines. We see that Eqs. (26) and (27) yield a better fit of the numerical results for large values of the longitudinal density ratio or for small values of the average density, which are the expected ranges of values of prominence threads.

Taking into account the results displayed in Fig. 8, the ratio τ_D, 0/τ_D, 1 could be used as another seismology tool to compute the parameters χ or ⟨ρ_i⟩/ρ_i, 0 of threads and provide results that are independent of the ones obtained through the use of the period ratio. Nevertheless, the true efficiency of this ratio and the approximations given by Eqs. (26) and (27) should be explored in a future work. Soler et al. (2015) and Ruderman et al. (2016) demonstrated that the fundamental variable that determines the value of P₀/P₁ is the ratio ⟨ρ_i⟩/ρ_i, 0, with no dependence from the specific longitudinal density profile of the thread. A similar general behaviour has not yet been established for the ratio τ_D, 0/τ_D, 1, which might depend on other additional factors that would constrain the range of applicability of the approximations presented here.

3.2.5. Amplitude of the oscillations

Figure 9 shows how the average density of the tube affects the amplitudes of the transverse oscillations of the threads (left panel) and their periods (right panel).

Fig. 9. Displacement (left) and oscillation period (right) as functions of average density. Black crosses correspond to simulations of tubes with uniform density and an amplitude of the initial perturbation given by V0 = 0.05cA. Red symbols represent the results with a Lorentzian profile and V0 = 0.05cA. Blue diamonds represent simulations of tubes with uniform density but an initial kinetic energy that is the same as the cases with the Lorentzian profile.

Firstly, we compare two series of simulations in which the amplitude of the initial perturbation is V₀ = 0.05c_A, in and the fundamental mode of oscillation, n = 1, is considered: black crosses correspond to uniform tubes with the same density as the average density of the Lorentzian profiles, while red asterisks show the results for tubes with Lorentzian profiles. For a given average density, a tube with uniform density oscillates with a larger amplitude and shorter period (as already shown by Fig. 7). In addition, we see that for the case of the Lorentzian profile the amplitude of the oscillation increases with the average density, while for a uniform thread it remains almost constant.

These results for the longitudinally homogeneous tube are in good agreement with the linear analysis for the transverse kink mode in the thin tube approximation performed by Ruderman & Goossens (2014) and Terradas et al. (2018), who showed that the displacement of the tube is given by

$$\begin{array}{c}\hfill A_0=\frac{L_z{V}_0}{\pi R{c}_{\mathrm{A},\mathrm{in}}}\sqrt{\frac{1+{\rho }_{\mathrm{ex}}/{\rho }_{\mathrm{i},0}}{2}}=\frac{L_z\alpha }{\pi R}\sqrt{\frac{1+{\rho }_{\mathrm{ex}}/{\rho }_{\mathrm{i}}}{2}}.\end{array}$$(28)

In the above formula, we used the fact that in our simulations the amplitude of the initial perturbation is given by V₀ = αc_A, in. Therefore, the displacement is almost a constant since ρ_i, 0 ≫ ρ_ex and the rest of parameters have been kept fixed for these series of simulations. The physical reason behind this behaviour is that all the simulations of this series start with the same total kinetic energy, independently from the average density of the tube.

We note that in a strict sense, Eq. (28) is only applicable to the linear regime where A₀ ≪ 1, while our simulations present much larger values of the displacement. However, it seems that the predicted trend also appears in the non-linear regime (at least for the set of parameters explored in the present work).

The behaviour of the simulations with a Lorentzian profile can be explained as follows. Although the value of χ varies throughout the series, the density at the central part of the tube and, consequently, the internal Alfvén speed remain fixed, which implies that the amplitude of the initial velocity perturbation is also fixed. Therefore, as the average density of the tube is increased (or χ is decreased), the total kinetic energy rises, which produces oscillations of larger amplitude.

In addition, for a given value of the averaged density, the internal Alfvén speed of the tube with a uniform density is larger than the one with a Lorentzian profile, since ρ_i, 0 is smaller in the former case. This is the reason why all the simulations with uniform density show larger displacements.

Next, we compared simulations with different density profiles but with the same initial kinetic energy. For that, we performed a new series of simulations of threads with uniform density but varied the amplitude of the initial perturbation, V₀, so their total kinetic energy would match that of the threads with a Lorentzian profile. The results of this new series are shown in Fig. 9 as blue diamonds. Surprisingly, in the left panel we see that the results of both series still differ, with the amplitudes for the cases with Lorentzian profiles being slightly larger.

An explanation for the mentioned discrepancy could be that, although the total kinetic energy is the same, a larger fraction of it is present in the central region of the tube in the case of the Lorentzian profile in comparison with the case of the uniform profile. To check the validity of this explanation, we computed the amplitude of the oscillation at a different position of the tube, namely z = −0.25L_z. At that position, for a given value of ⟨ρ_i⟩/ρ_i, 0, the uniform profile has a larger density (and larger kinetic energy) than the Lorentzian profile. Therefore, smaller amplitudes would be expected for the latter. The results of these computations are shown in Fig. 10. It can be seen that the amplitudes for the cases with a uniform density are slightly larger, which agrees with the explanation provided.

Fig. 10. Same as Fig. 9, but for z = −Lz/4.

Finally, the right panels of Figs. 9 and 10 show that the results of the two series of simulations of threads with homogeneous density overlap. Thus, we find that in these simulations the periods of the oscillations are not affected by the amplitude of the initial perturbation, in agreement with the behaviour predicted by Eq. (20). However, we remind the reader this formula is only strictly applicable to the linear regime, that is, for oscillations with amplitudes A₀ ≪ 1. Figures 9 and 10 show that the simulated oscillations have values of A₀ that go from ∼0.15 to ∼0.5, which means that they are out of the range of applicability of Eq. (20). According to Ruderman & Goossens (2014), non-linearity produces an increase of the oscillation frequency proportional to the square of the amplitude A₀. More precisely, the relation between the frequencies in the non-linear and linear regimes is given by ${\omega }_{\mathrm{nonl}}={\omega }_{\mathrm{lin}}(1+A_0^2\mathrm{{\rm Y}})$, where Υ is a small parameter that depends on the ratio ρ_i/ρ_ex and the geometry of the flux tube. Therefore, smaller periods are obtained as the amplitude is increased. Nevertheless, Ruderman & Goossens (2014) also found that the non-linear frequency shift is usually small and it would become noticeable after a large number of periods. Hence, even larger amplitudes than the ones used here would be required to obtain a strong variation in the oscillation periods.

4. Forward modelling of inhomogeneous flux tubes

We performed forward modelling of our simulations to compute the synthetic intensity of the Hα line emerging from the oscillating thread. We chose the Hα spectral line for the present study because it is associated with the neutral component of the plasma (and, consequently, related to low temperatures), and it is typically used for imaging the cool prominence fine structures (Lin 2011).

To obtain the synthetic emission of the Hα line, we implemented the approximate method developed by Heinzel et al. (2015). We used the reference values for a prominence at a height of 10000 km provided by Table 1 of Heinzel et al. (2015), so the plasma in the interior of the thread has a temperature of T = 10 000 K, a pressure of P = 0.005 Pa, and an ionisation degree of i = 0.7. The lighter external plasma is fully ionised, i = 1, and initially has the same pressure as the internal plasma in order to fulfill the condition of mechanical equilibrium. Since the internal and external densities follow the relation ρ_i, 0 = 100ρ_ex, the initial temperature of the external medium is T ≈ 8.5 × 10⁵ K. As time advances, at every step of the simulation the ionisation degree of each point of the numerical domain is computed by performing a bilinear interpolation from the data included in Table 1 of Heinzel et al. (2015).

In this section, we increase the resolution of our simulations and use a numerical domain of 400 × 400 × 100 points instead of the 200 × 200 × 100 described in the previous section. Increasing the numerical resolution provides a more accurate description of the non-linear dynamics of the flux tube. It allows us to better resolve smaller scales that were not captured by the low-resolution simulations, as we show in more detail in the paragraphs below, and obtain better predictions of the results of the observations.

Before presenting the Hα synthetic profiles, for reference purposes we show density snapshots taken from simulations with V₀ = 0.05c_A and n = 1 in Fig. 11. The snapshots correspond to horizontal cuts at z = 0, where the displacement is at its maximum. The top row shows the results for the simulation with χ = 10, while the bottom row corresponds to χ = 100. We see that up to the time t/t₀ = 400, the tube keeps its cylindrical shape almost unaltered. On the contrary, at later times its boundary is strongly distorted as the KHI vortexes develop. For detailed investigations on the onset and evolution of the KHI in magnetic flux tubes, we refer the reader to the works of, for example, Terradas et al. (2008a, 2018), Antolin et al. (2014, 2017), Magyar et al. (2015), Magyar & Van Doorsselaere (2016), or Díaz-Suárez & Soler (2021).

Fig. 11. Density colour maps at z = 0 from simulations with V0 = 0.05cA, n = 1, and a numerical domain of 400 × 400 × 100 points. Top panel: case with χ = 10; bottom panels: case with χ = 100. Dashed lines represent the density contour ρ = 50ρref at t = 0. An animation of this figure is available online.

The top panels of Fig. 11 can be compared with those of Fig. 3. They represent oscillations of threads with the same physical parameters (density profile, longitudinal oscillation mode, and amplitude of the perturbation) but differ in the numerical resolution employed in each simulation; the simulation depicted in Fig. 3 used a lower numerical resolution, while the simulation represented in Fig. 11 used a higher resolution. This comparison shows the importance of using a higher numerical resolution (a smaller cell size) to properly simulate the evolution of the KHI. It allows us to resolve smaller scales, associated with larger azimuthal wave numbers, which have larger growth rates (Soler et al. 2010b), and show more details of the development of the fine structure during the thread oscillation, as demonstrated by Antolin et al. (2015). These authors also found that improving the numerical resolution produces an increase of small-scale vortexes without altering the global oscillation of the tube.

Then, in Figs. 12 and 13 we present the results of applying the approximate forward modelling to our simulations. We took into account the whole 3D simulation to compute the Hα synthetic profile and show how the full longitudinal structure of the flux tube would be seen in the observations, revealing in that way the effect of the inhomogeneity in density. We considered two different lines of sight (LOS): the first one, denoted as x_LOS, is parallel to the direction of oscillation, which goes along the x-axis; the second one, y_LOS, is perpendicular to the direction of oscillation, that is, along the y-axis. Both LOS are perpendicular to the longitudinal axis of the tube. In each of the panels of Figs. 12 and 13 the top and bottom panels correspond to the x_LOS and the y_LOS, respectively. In addition, we have added density contours for the values ρ = 1.05ρ_ex and ρ = 50ρ_ex as white dashed lines. The density contours are computed at the plane that is perpendicular to the LOS and crosses the longitudinal axis of the tube, that is, the plane x = 0 for the x_LOS and the plane y = 0 for the y_LOS. These lines serve to illustrate how the whole magnetic tube, and not only the dense central part, oscillates.

Fig. 12. Intensity in Hα line at different times of the oscillation for the case with χ = 10. In each snapshot, the top panel corresponds to the xLOS, while the bottom panel corresponds to the yLOS. White dashed lines show density contours for ρ = 1.05ρex and ρ = 50ρex. An animation of this figure is available online.

Fig. 13. Same as Fig. 12, but for the simulation with χ = 100. An animation of this figure is available online.

We also applied the approximate forward modelling method to a simulation of the mode n = 2 and the longitudinal density ratio χ = 10. The rest of the parameters of the equilibrium state and the initial perturbation are the same as the ones used in the previous simulations. Figure 14 presents several snapshots from this simulation. The density contour ρ = 1.05ρ_ex approximately shows how the whole magnetic flux tube oscillates, with its maximum displacement occurring at z = −12.5L₀ and z = 12.5L₀.

Fig. 14. Intensity in Hα line at different times of the oscillation for the case with χ = 10 and n = 2 (first harmonic). In each snapshot, the top panel corresponds to the xLOS, while the bottom panel corresponds to the yLOS. White dashed lines show density contours for ρ = 1.05ρex and ρ = 50ρex. An animation of this figure is available online.

The comparison between Figs. 12 and 14 shows how important the effect of the longitudinal distribution of density in the thread is from the observational point of view. The central part of the tube appears much brighter than the remaining regions. The reason is that the intensity has an approximately quadratic dependence on the density, so the variations of density are very strongly reflected on the intensity profiles. It also reveals one of the main practical problems that the field of prominence seismology has to face when using tools like the period ratio P₀/P₁. On one hand, the Hα line provides information about the cool and dense central part of the flux tube and it allows to measure the fundamental mode of oscillation. On the other hand, the first harmonic has a node at z = 0 while its maximum displacement takes place in the dark regions of the Hα profile (see Fig. 14). Therefore, another filter associated with the lighter and hotter plasma would be needed to measure the properties of this oscillation mode.

Another remarkable feature of the intensity profiles is the development of bright and dark fringes as time advances. These fringes are the observational signatures of the deformation of the tube boundary caused by the KHI. The stripes are more evident when the thread is observed along the LOS that is perpendicular to the direction of oscillation. However, for the case with χ = 10, they can also be seen in the other LOS. For χ = 100, there are variations of brightness along the x_LOS, but the stripes are not so easily discerned. The development of this fine strand-like structure as a consequence of the KHI was already described by Antolin et al. (2014, 2016, 2017) for the case of coronal loops and by Antolin et al. (2015) for the case of prominences.

There are two main differences between the results we present in this work and the results shown by Antolin et al. (2014, 2015, 2016, 2017). Firstly, these authors performed forward modelling of spectral lines associated with higher temperatures, while we used a line related to cool plasma. Therefore, our forward modelling results and theirs can be compared qualitatively but not quantitatively. The other difference is that Antolin et al. (2014, 2015, 2016, 2017) considered a tube model with no longitudinal variation of density. They found that the KH vortexes are present at every longitudinal position of the tube. We expect the same to occur in the case with the longitudinal inhomogeneity, although they can only be seen in the intensity profiles near z = 0, where the plasma is denser. By inspecting density snapshots at different positions of the flux tube (which we do not present here), we checked that this is indeed what occurs.

In the last place, Figs. 15 and 16 show how a horizontal cut at z = 0 would be seen as a function of time for the cases with χ = 10 and χ = 100, respectively. Again, the parallel and perpendicular lines of sight are displayed in the top and bottom panels, respectively. The oscillatory motion of the thread is more clearly appreciated in the perpendicular LOS, y_LOS, which allows us to easily trace the displacement of the central part of the tube, the development of the instability, and the attenuation of the oscillation due to the resonant absorption. From these panels, we can also estimate that the width of the most conspicuous bright and dark strands is of the order of w ∼ R/5 or R/4 (that is, from 200 to 250 km), which corresponds to an angular diameter of δ ∼ 0.3″. Therefore, these fine structures could be resolved by Hα instruments such as the CRisp Imaging SpectroPolarimeter (CRISP Scharmer et al. 2008) installed at the Swedish Solar Telescope (SST, Scharmer et al. 2003), or the Visible Imaging Spectrometer (VIS, Cao et al. 2010) installed at the Goode Solar Telescope (GST, Goode & Cao 2012), which have spatial resolutions of the order of ∼ 0.1″(∼ 70 km) (see e.g., Jing et al. 2019; Froment et al. 2020). However, we note that, as already mentioned in Sect. 2.1, the value of R used in this work is slightly larger than what has been observed for prominence threads. Typically, the radius of a thread ranges from 50 km to 500 km (Okamoto et al. 2007; Lin et al. 2008; Lin 2011), which would give a range of 10–125 km for the bright and dark fringes. Consequently, those instruments would not be able to resolve the KH-strands of thin threads, but they would allow us to observe the fine structure of thicker threads. More recent instruments, such as the Visible Broadband Imager (VBI, Wöger et al. 2021) installed at the Daniel K. Inouye Solar Telescope (DKIST, Rimmele et al. 2020), should be able to detect strands with a width of the order of ∼ 25 km.

Fig. 15. Intensity of Hα line at z = 0 as a function of time from a simulation with a density ratio of χ = 10.

Fig. 16. Intensity of Hα line at z = 0 as a function of time from a simulation with a density ratio of χ = 100.

On the other hand, when observed along the direction of oscillation, the amplitude of this motion cannot be directly determined from these time-distance diagrams. Nonetheless, some evidence of the oscillation can still be found. Apart from the bright and dark fringes that can be associated with the KHI, there are also periodic variations of brightness. These variations are related to the contractions and expansions that appear in the tube before the KH vortexes develop. As shown by Antolin et al. (2017) and Terradas et al. (2018), the width of the tube periodically increases in the x-direction and decreases in the y-direction, and vice versa. A contraction in one direction produces a decrease in the Hα intensity along the parallel LOS but an increase in the perpendicular LOS, while expansions have the opposite effect. This behaviour can be checked by comparing the top panels of Figs. 15 and 16 with their corresponding bottom panels.

By comparing Figs. 15 and 16, it can be seen that the thread with a lower density ratio oscillates with a larger amplitude and a longer period, in agreement with the results presented in Sect. 3.2.1. In addition, it presents larger variations of brightness in its x_LOS profile, which is a consequence of having a larger amplitude of the oscillation. In both cases, the fringes associated with the non-linear stage of the KHI start to be clearly seen at around the same time step, t/t₀ ≈ 550, which is after 1 oscillation period for χ = 10 and about 1.5 periods for χ = 100. It is also interesting to note how the density contours follow the oscillation of the bright region of the tube during the first periods, but they become out of phase as time advances. This is a consequence of the deformations caused by the KHI and the way the contours and the Hα intensity are computed: the former uses the values of density at only one plane, while the latter is the result of integrating along the whole LOS (which crosses different planes with different values of density).

5. Summary and future work

In the present work, we performed a numerical study of the longitudinal fundamental and first harmonic modes of transverse oscillations in inhomogeneous threads. We confirmed the findings of Soler et al. (2015) showing that, if the central density of the tube is kept fixed, the oscillation period is reduced as the ratio between the densities at the centre of the tube and at its ends, χ, is increased. We also confirm that the period ratio of the two modes, P₀/P₁, increases as the longitudinal density ratio is increased, or as the average density of the tube is decreased.

Then, we analysed the damping times of the oscillations. The attenuation is mainly due to the process of resonant absorption (see e.g., Ionson 1978; Ruderman & Roberts 2002). There is also a small contribution due to the non-linear generation of modes with a high azimuthal wave number (Ruderman & Goossens 2014; Soler 2017) and the development of the KHI (Terradas et al. 2018; Van Doorsselaere et al. 2021). We find that the damping times also decrease as the longitudinal density ratio of the tube is increased. However, the damping-time-to-period ratio increases with χ, which means that the increase of χ causes a relatively larger reduction of the oscillation periods in comparison with the damping times. This result qualitatively agrees with those of Arregui et al. (2011), who considered a longitudinal density structuring and found a slight increase of the damping ratio for very short threads (which in our model would be represented by an increasing value of χ).

In addition, a more pronounced effect of the increase of χ on the damping ratio is found for the case of the first harmonic. The reason is that the efficiency of resonant absorption depends on the density contrast between the internal and external densities of the flux tube and the thickness of the radial transition layer (Ruderman & Roberts 2002; Arregui et al. 2008; Soler et al. 2009). In the present work, the density contrast varies with height and with the longitudinal density ratio χ. At the height of the maximum displacement of the first harmonic, there are larger variations of the density contrast with χ than at the height where the fundamental mode has its maximum displacement.

Taking into account the results described in the previous paragraphs, we proposed using the ratio τ_D, 0/τ_D, 1 in the same way the period ratio P₀/P₁ is used for prominence seismology. From the solutions of the 2D eigenvalue problem, we obtained two approximate expressions, Eqs. (26) and (27), for the dependence of τ_D, 0/τ_D, 1 on the parameters χ and ⟨ρ_i⟩/ρ_i, 0, respectively. However, the utility of this parameter for extracting information from the observations should be analysed in more detail in a future work. For instance, in contrast with the period ratio, which only depends on the relation between the average density and the density at the central part of the tube, as shown by Soler et al. (2015) and Ruderman et al. (2016), the ratio τ_D, 0/τ_D, 1 might depend on the specific longitudinal profile of density and its application might not be as straightforward. Nevertheless, we expect that this longitudinal dependence would be less relevant than its dependence on the profile of the radial transition layer.

The last part of the present work consisted of performing forward modelling from our 3D simulations. We computed synthetic Hα profiles using the approximate method developed by Heinzel et al. (2015). In this way, we showed the effect that the Lorentzian profile of density would have on the observations: only the central denser part of the thread would be visible in the intensity profiles, while most of the tube would appear dark in comparison. This is an important issue for the field of prominence seismology. The Hα profiles allow us to detect and measure the characteristics of the fundamental oscillation mode with relative ease, since its maximum displacement happens at the bright part of the thread. On the contrary, the first longitudinal harmonic has a node in that bright central region, and the maximum displacement occurs in the dark regions, making it extremely difficult to measure the properties of this mode. We recall that, to the best of our knowledge, no unequivocal detection of the first harmonic in prominence threads has been reported to date. This issue can be resolved by combining the Hα observations with measures in spectral lines associated with the emission from lighter and hotter plasma.

Apart from the periods and damping times of the oscillations and the density distribution of the thread, the intensity profiles also present some observational signatures of non-linear effects that take place as the tube oscillates. The large amplitudes of the initial perturbations used in our simulations cause the appearance of strong velocity shears that trigger the KHI (Browning & Priest 1984; Terradas et al. 2008a). This instability generates large deformations of the external layers of the tube, as shown by Antolin et al. (2014, 2017), Magyar et al. (2015), Magyar & Van Doorsselaere (2016), or Terradas et al. (2018). These deformations are reflected in the Hα profiles by longitudinal bright and dark stripes. This is the same mechanism that can generate the strand-like structure seen in coronal loops (Antolin et al. 2014, 2016, 2017) and prominences (Antolin et al. 2015). We estimate that, for a typical prominence thread, the size of this fine structure varies in the range from 10 to 125 km. Therefore, it should be possible to observe these fine strands in thick threads using instruments such as SST/CRISP (Scharmer et al. 2003, 2008) and GST/VIS (Cao et al. 2010; Goode & Cao 2012), which have an Hα spatial resolution of ∼ 70 km, or DKIST/VBI (Rimmele et al. 2020; Wöger et al. 2021) (with a resolution of ∼ 25 km).

Here, we describe the dynamics of magnetic flux tubes through the ideal MHD equations. However, prominence threads are made of partially ionised plasma. Thus, for a more realistic analysis, the effects of the interaction between the ionised and neutral components of the plasma should be taken into account. Soler et al. (2009, 2014) showed that ambipolar diffusion has a negligible impact on the attenuation of transverse oscillations in comparison with the mechanism of resonant absorption. Nevertheless, the ion-neutral interaction may have an important effect on non-linear processes, such as the heating of the plasma, specially as the KHI develops and smaller and smaller scales are generated, in which the energy dissipation due to ion-neutral collisons is more efficient (Leake et al. 2005; Zaqarashvili et al. 2011; Soler et al. 2013; Popescu Braileanu et al. 2019). In addition, this energy dissipation may oppose the action of the ponderomotive force, as shown by Martínez-Gómez et al. (2018) and Ballester et al. (2020) for the case of non-linear Alfvén waves.

Other effects that influence the oscillations in flux tubes but that we did not consider are the presence of mass flows, which have already been studied using uniform or piecewise constant models (Terradas et al. 2008b; Soler & Goossens 2011; Erdélyi et al. 2014; Bahari et al. 2020; Sadeghi et al. 2021) but not with longitudinally inhomogeneous density profiles, or the twist, curvature, and cross-section of the tube, which have been analysed for coronal loops (Verth & Erdélyi 2008; Verth et al. 2008; Ruderman et al. 2008; Morton & Erdélyi 2009; Karami & Bahari 2012) but not for prominence threads. In addition, Martínez-Gómez et al. (2015) showed that since the plasma is partially ionised in prominence threads, they are always unstable in the presence of longitudinal shear flows, which trigger the KHI. The inclusion of some of the aforementioned effects in our model is left for future works.

Movies

Movie 1 associated with Fig. 11 (Fig11_top) Access here

Movie 2 associated with Fig. 11 (Fig11_bottom) Access here

Movie 3 associated with Fig. 12 Access here

Movie 4 associated with Fig. 13 Access here

Movie 5 associated with Fig. 14 Access here

Acknowledgments

D. M. and E. K. acknowledge support from the European Research Council through the Consolidator Grant ERC-2017-CoG-771310-PI2FA. D. M. also acknowledges support from the Spanish Ministry of Science and Innovation through the grant CEX2019-000920-S of the Severo Ochoa Program. This publication is part of the R+D+i project PID2020-112791GB-I00, financed by MCIN/AEI/10.13039/501100011033. We also thank the anonymous referee for constructive remarks and helpful suggestions.

References

All Figures

	Fig. 1. Sketch of equilibrium state for a thread with a density ratio of χ = 10. Density contours for ρ = 25ρex (red), ρ = 50ρex (green), ρ = 75ρex (dark blue), and ρ = 100ρex (light blue). Magnetic field lines, with line-tying conditions applied at the top and bottom planes, are shown in grey.
In the text

	Fig. 2. Density profiles of equilibrium state. Top panel: vertical profiles for several values of the density ratio: χ = 5 (black), χ = 10 (green), χ = 20 (light blue), χ = 50 (dark blue), and χ = 100 (red). Bottom panel: radial profile at z = 0.
In the text

	Fig. 3. Snapshots of a simulation with χ = 10, V0 = 0.05cA, in, n = 1, and a numerical domain of 200 × 200 × 100 points. Top: density colour maps at z = 0. Bottom: density colour maps at y = 0. Dashed lines represent the density contour ρ = 50ρref at t = 0.
In the text

	Fig. 4. Snapshots of a simulation with χ = 10 and n = 2. Top: density colour maps at z = −Lz/4 = −12.5L0. Bottom: density colour maps at y = 0. Dashed lines represent the density contours ρ ≈ 15ρref (top panels) and ρ = 50ρref (bottom) at t = 0.
In the text

	Fig. 5. Displacement of tube axis at different longitudinal positions from simulations of modes n = 1 (top) and n = 2 (bottom). The solid and dotted lines represent oscillations with periods given by Eq. (20) for tubes of uniform density equal to ρi, 0 and ⟨ρi⟩, respectively. The dashed lines represent the exponential decay of the displacement of the fundamental mode.
In the text

Fig. 6. Results of the parametric study of transverse oscillations. Top left: displacement of centre of mass as a function of time from simulations with n = 1 and with χ = 5 (black crosses), χ = 10 (green asterisks), χ = 20 (light blue diamonds), χ = 50 (dark blue triangles), and χ = 100 (red squares). Top right: oscillation period as a function of the density contrast, χ. Bottom left: damping time as a function of χ. Bottom right: ratio of damping time to period as a function of χ. In the top left panel, solid lines represent fits of the numerical results. In the remaining panels, the black solid and dashed lines correspond to the results obtained by solving the 2D eigenvalue problem for n = 1 and n = 2, respectively. Symbols represent the results from the simulations. For n = 1, Pref = P0; for n = 2, Pref = P1.

In the text

	Fig. 7. Period ratio as a function of density ratio (left) and as a function of the average density of the tube (right). Dashed lines correspond to the case of tubes with uniform density. Dotted lines represent the approximation given by Eq. (24).
In the text

Fig. 8. Damping time ratio as a function of density ratio (left) and as a function of the average density of the tube (right). Red asterisks represent the results from simulations. Solid black lines correspond to the results from the eigenvalue problem. Horizontal dashed lines represent the case of tubes with uniform density. The red dotted lines represent the approximations given by Eqs. (24) and (25). The blue dotted-dashed lines in the left and right panels represent the fits given by Eqs. (26) and (27), respectively.

In the text

Fig. 9. Displacement (left) and oscillation period (right) as functions of average density. Black crosses correspond to simulations of tubes with uniform density and an amplitude of the initial perturbation given by V0 = 0.05cA. Red symbols represent the results with a Lorentzian profile and V0 = 0.05cA. Blue diamonds represent simulations of tubes with uniform density but an initial kinetic energy that is the same as the cases with the Lorentzian profile.

In the text

	Fig. 10. Same as Fig. 9, but for z = −Lz/4.
In the text

	Fig. 11. Density colour maps at z = 0 from simulations with V0 = 0.05cA, n = 1, and a numerical domain of 400 × 400 × 100 points. Top panel: case with χ = 10; bottom panels: case with χ = 100. Dashed lines represent the density contour ρ = 50ρref at t = 0. An animation of this figure is available online.
In the text

	Fig. 12. Intensity in Hα line at different times of the oscillation for the case with χ = 10. In each snapshot, the top panel corresponds to the xLOS, while the bottom panel corresponds to the yLOS. White dashed lines show density contours for ρ = 1.05ρex and ρ = 50ρex. An animation of this figure is available online.
In the text

	Fig. 13. Same as Fig. 12, but for the simulation with χ = 100. An animation of this figure is available online.
In the text

	Fig. 14. Intensity in Hα line at different times of the oscillation for the case with χ = 10 and n = 2 (first harmonic). In each snapshot, the top panel corresponds to the xLOS, while the bottom panel corresponds to the yLOS. White dashed lines show density contours for ρ = 1.05ρex and ρ = 50ρex. An animation of this figure is available online.
In the text

	Fig. 15. Intensity of Hα line at z = 0 as a function of time from a simulation with a density ratio of χ = 10.
In the text

	Fig. 16. Intensity of Hα line at z = 0 as a function of time from a simulation with a density ratio of χ = 100.
In the text