© ESO 2021

1. Introduction

The solar atmosphere shows a variety of dynamic structures that are much denser than the background corona. These formations, guided by the curved magnetic field and most likely created by thermal instabilities, have their own dynamics and are undoubtedly affected by the gravity force. Examples of such configurations are solar prominences or coronal rain among others.

A question that immediately arises regarding these structures, and especially for solar prominences, is how they can be suspended against gravity above the photosphere for long periods. This has led to the development of theoretical prominence models that include magnetic dips where the gravity force projected along the field lines is zero. Analytical examples of such equilibrium configurations were developed some years ago by Hood & Anzer (1990), Fiedler & Hood (1992), and de Bruyne & Hood (1993). Later, equilibrium models in flux rope configurations were obtained by Low & Zhang (2004), Petrie et al. (2007), and Blokland & Keppens (2011). More recently, numerical equilibria are achieved through relaxation processes by solving the time dependent magnetohydrodynamics (MHD) equations in two-dimensional (2D) configurations (Hillier & van Ballegooijen 2013; Terradas et al. 2013; Luna et al. 2016a). Three-dimensional (3D) stable magnetic configurations with dips are numerically investigated by Terradas et al. (2015, 2016) and Adrover-González & Terradas (2020). In these last works the presence of dips is not crucial for the stability of the prominence. However, it seems that there is the common idea in the solar community that threads or blobs automatically fall along the field lines because of gravity force and that if threads are suspended it is because of magnetic dips. However, from the force-balance point of view, the gas pressure might still support a short segment of prominence as numerically illustrated by Kohutova & Verwichte (2017a) and Mackay & Galsgaard (2001).

A characteristic of the Hood & Anzer (1990) model is that the equilibrium is in general unstable to lateral displacements, meaning that the dense and cold part of the prominence body tends to fall along magnetic field lines essentially owing to the gravity force (see also the numerical experiments of Kohutova & Verwichte 2017a). This was demonstrated in the stability analysis of de Bruyne & Hood (1993). Therefore, the presence of magnetic dips in the configuration is not a warranty for stability. On the contrary, in the magnetic configuration studied numerically in Terradas et al. (2013) the gravitational instability reported in the Hood & Anzer (1990) model was clearly not taking place. This is the main motivation of part I of the present work, namely to understand the reasons for gravitational instability in a global sense, and not necessarily under the presence of magnetic dips. The effect of magnetic dips will be addressed in part II.

At first glance, this instability may seem to be directly related to magnetic buoyancy, which is a phenomenon that has been investigated in great detail by Parker (1979). But it turns out that magnetic buoyancy is not related to this kind of process since it also appears in the absence of fluctuations in the magnetic field so the instability is purely gravitational.

In this work we address gravitational instability using a very simple configuration with the consideration of gravity. This allows us to disentangle the main physical processes that occur in the system. We consider a circular magnetic flux tube in which gravity has a spatially varying parallel component along the field lines. We use elementary physics describing a plasma in this configuration, that is the interplay between gas pressure and gravity force. The magnetic field is assumed not to change as a result of the presence of the density enhancement or thread and simply acts as a guide for longitudinal motions. We do not include non-ideal effects, which are most likely important in the solar corona (and chromosphere), because in our opinion the purely mechanical situation under line-tying conditions in the absence of dissipative mechanisms has received little attention in the literature and needs to be understood more deeply.

2. Background equilibrium solution

We start with the description of our model. We assume that the shape of the magnetic field is circular and that its cross-section, A₀, is constant along the field (see dotted curves in Fig. 1 top panel). We denote by s the coordinate along the tube, starting from the left foot. The total length of the magnetic tube is L and the top or apex is located at s = L/2. Gravity is pointing downwards and the projection along the field lines in this circular configuration is given by

Fig. 1. Sketch of the basic model. Upper panel: situation in which the plasma thread, represented by the dark rectangle of length lt, is located at the tube apex. This configuration is in an equilibrium that can be stable or unstable. Lower panel: new equilibrium of the thread at a lower height. The curved magnetic field is uniform and directed along the tube axis. The gravity force is pointing downwards and is projected along the field line according to Eq. (1). The footpoints are supposed to be anchored at the base of the corona, represented by the shaded grey area.

$$\begin{array}{c}\hfill g_{\Vert }(s)=-gcos\left(\pi \frac{s}{L}\right),\end{array}$$(1)

where the radius of the field line is R = L/π.

The hydrostatic equation along the tube is simply

$$\begin{array}{c}\hfill \frac{\mathrm{d}p}{\mathrm{d}s}(s)=\rho (s)g_{\Vert }(s),\end{array}$$(2)

meaning that the pressure derivative with position must balance the projected gravity force along the field lines. Since we assume an isothermal atmosphere the sound speed, c_s0, is constant and we write

$$\begin{array}{c}\hfill p(s)=\frac{c_{\mathrm{s}0}^2}{\gamma }\rho (s),\end{array}$$(3)

where γ = 5/3. Introducing this expression into Eq. (2) we find the following simple differential equation

$$\begin{array}{c}\hfill \frac{\mathrm{d}\rho }{\rho }=\frac{\gamma }{c_{\mathrm{s}0}^2}{g}_{\Vert }(s)\mathrm{d}s,\end{array}$$(4)

which is integrated directly to give

$$\begin{array}{c}\hfill \rho (s)={\overline{\rho }}_{0}exp[-\frac{\gamma g}{c_{\mathrm{s}0}^2}\frac{L}{\pi }sin\left(\pi \frac{s}{L}\right)],\end{array}$$(5)

where ${\overline{\rho }}_0$ is the reference density at s = 0 and s = L, that is at the footpoints of the magnetic tube. According to Eq. (3) gas pressure and density have the same dependence with s.

3. Equilibrium solutions containing threads

A stratified prominence unbounded with height, as for example in the Hood & Anzer (1990) model, can be seen as a continuous collection of threads stacked in the vertical direction but located at magnetic dips. It is crucial to understand the stability of a density enhancement or equivalently a thread in the solar atmosphere, even in the absence of magnetic dips as we explained in Sect. 1.

The idea is to analyse the simple situation of a finite mass (called a thread, hereafter) with density ρ_t located on the curved magnetic field (see Fig. 1 for a sketch of the model). The question is whether this equilibrium is stable or unstable and we investigate the conditions in which the system allows equilibrium solutions at lower heights.

The basic equation of our system is the momentum equation

$$\begin{array}{c}\hfill {\rho }_{\mathrm{t}}\frac{{\mathrm{d}}^{2}s}{\mathrm{d}{t}^2}=-\frac{\mathrm{d}p}{\mathrm{d}s}(s)+{\rho }_{\mathrm{t}}{g}_{\Vert }(s),\end{array}$$(6)

where s is the coordinate of the centre of the thread along the curved path. Since the thread is surrounded by plasma, fluctuations in gas pressure at the edges of the thread in the longitudinal direction create a net force. We need to provide an expression for the gas pressure derivative. The other force existing in the system is the gravity projection along the magnetic field line given in Eq. (1). Before we go into the details about the calculation of the gas pressure derivative, it is necessary to describe the general features of our system from a mathematical point of view.

In order to analyse the stability of the system, we write Eq. (6) as the vector field

$$\left(\begin{array}{c}\dot{s}\\ \dot{v}\end{array}\right)=\left(\begin{array}{c}v\\ -\frac{1}{{\rho }_{\text{t}}}\frac{\text{d}p}{\text{d}s}(s)+g_{\parallel }(s)\end{array}\right),$$(7)

where the dot represents the time derivative. The system has equilibrium or fixed points when

$$\begin{array}{cc}& v=0,\hfill \\ \hfill & -\frac{1}{{\rho }_{\mathrm{t}}}\frac{\mathrm{d}p}{\mathrm{d}s}(s)+g_{\Vert }(s)=0,\hfill \end{array}$$(8)

and the corresponding solution for the position is denoted as s_e; the velocity at the equilibrium point is always zero. The linearised version of Eq. (7) is written as

$$\left(\begin{array}{c}\dot{s}\\ \dot{v}\end{array}\right)=A\left(\begin{array}{c}s\\ v\end{array}\right),$$(9)

where the matrix A is the Jacobian; and

$$\begin{array}{c}\hfill \mathbf{A}=\left(\begin{array}{cc}\frac{\partial \dot{s}}{\partial s}& \frac{\partial \dot{s}}{\partial v}\\ \frac{\partial \dot{v}}{\partial s}& \frac{\partial \dot{v}}{\partial v}\\ \end{array}\right)=\left(\begin{array}{cc}0& 1\\ f(s)& 0\\ \end{array}\right),\end{array}$$(10)

where

$$\begin{array}{c}\hfill f(s)=\frac{\mathrm{d}}{\mathrm{d}s}[-\frac{1}{{\rho }_{\mathrm{t}}}\frac{\mathrm{d}p}{\mathrm{d}s}(s)+g_{\Vert }(s)]·\end{array}$$(11)

The eigenvalues of A at the equilibrium points (s = s_e) provide valuable information about their stability. If f(s_e) > 0, we have

$$\begin{array}{c}\hfill {\lambda }_{\pm }=\pm \sqrt{f(s_{\mathrm{e}})}.\end{array}$$(12)

In this case the fixed point is a saddle point because we have real values with opposite signs (e.g., Jordan & Smith 1987). This point is unstable.

If f(s_e) < 0, where the linearised system is exactly the same as that of the simple harmonic oscillator equation, the eigenvalues are written as

$$\begin{array}{c}\hfill {\lambda }_{\pm }=\pm i\sqrt{|f(s_{\mathrm{e}})|}.\end{array}$$(13)

This corresponds to a linear centre that is a stable point. The last case is when f(s_e) = 0; this is a degenerate case but still provides information about the stability of the equilibrium points.

3.1. Thread initially at the tube apex

We first analyse the situation of the thread initially located at the top of a curved magnetic field (see Fig. 1 top panel). This mass is in equilibrium because the gravity force is zero at s = L/2 (see Eq. (1)) and the gradient of the background pressure is zero at the tube apex (see Eqs. (3) and (5)). The goal is to obtain an approximation for the pressure derivative in Eq. (8) once the thread has moved from the initial equilibrium position. In our model we apply line-tying conditions, meaning that pressure fluctuations are not lost along the field lines at the base of the corona. There is no energy leakage through the footpoints and mass flows along the magnetic tube injected from the photosphere are not allowed. Therefore we have perfect reflection at these points because of the big difference in density between the photosphere and corona. It is important to remark that this specific boundary condition allows us to have suspended threads at a given height even in the absence of magnetic dips and that the presence of a chromospheric layer (not included in our model) may affect how the pressure is balanced along the tube.

We denote the length of the thread as l_t and by a simple geometrical reasoning l_t/2 < s < L − l_t/2 because the thread edge cannot connect with the base of the corona. Under equilibrium at s = L/2 the length of the external part of tube that contains the thread is l₀ = L/2 − l_t/2 along each loop leg (see Fig. 1 top panel). We have the simple relation L = 2l₀ + l_t. Now we assume that the thread has moved to another location along the magnetic field (see Fig. 1 bottom panel). The distance on the left side is l₁ now (different from l₀), while the distance on the right side is l₂. We have that l₁ = s − l_t/2, and l₂ = L − l₁ − l_t = L − s − l_t/2. Quiescent prominences are typically 100 − 200 Mm long (L) and are about a 100-fold denser than the surrounding corona (ρ_t/ρ₀). However, it is known that the entire bodies of prominences consist in a collection of individual short threads of plasma located in long field lines so that only a fraction of each field line (l_t) is filled with dense plasma. In this study a wide range of thread parameters are utilised to exemplify a vast sort of prominence threads (see Mackay et al. 2010; Engvold 2015, for a detailed description of solar prominences).

Since the system is assumed to be dissipationless, we consider that the evolution during the motion of the thread is adiabatic. The behaviour of the plasma in the evacuated parts of the tube, that is where the density is low, is similar to that of a piston because of the line-tying conditions at the footpoints. The gas pressure can be written as the background, stratified from gravity, plus a term that depends on the position of the thread. The key idea is to apply the adiabatic law to each part of the tube, but to this term only, removing the contribution from the background,

$$\begin{array}{cc}& p_1{\left(l_1{A}_1\right)}^{\gamma }=p_0{\left(l_0{A}_0\right)}^{\gamma },\hfill \\ \hfill & p_2{\left(l_2{A}_2\right)}^{\gamma }=p_0{\left(l_0{A}_0\right)}^{\gamma },\hfill \end{array}$$(14)

where p₀ represents the reference pressure before the thread has moved, changing the length of the evacuated part from l₀ to l₁ and pressure from p₀ to p₁ for the left side. For the right side we have changes from l₀ to l₂ and from p₀ to p₂. Since we are considering a perfectly circular geometry we have that A₀ = A₁ = A₂ and the cross-sectional areas cancel in the previous expressions. Now from Eq. (14) we obtain that

$$\begin{array}{cc}& p_1(s)=p_0{\left(\frac{l_0}{l_1}\right)}^{\gamma }=p_0{\left(\frac{L/2-l_{\mathrm{t}}/2}{s-l_{\mathrm{t}}/2}\right)}^{\gamma },\hfill \\ \hfill & p_2(s)=p_0{\left(\frac{l_0}{l_2}\right)}^{\gamma }=p_0{\left(\frac{L/2-l_{\mathrm{t}}/2}{L-s-l_{\mathrm{t}}/2}\right)}^{\gamma }.\hfill \end{array}$$(15)

At this point we impose the condition that the plasma thread has to be in static and stationary equilibrium, namely, no flows are allowed and the solution does not change in time. These conditions imply that Eq. (7) reduces to Eq. (8). To simplify things we assume that the thread does not change its shape and keeps its length, l_t, constant. Now we approximate the gas pressure derivative that appears in Eq. (8) as

$$\begin{array}{c}\hfill \frac{\mathrm{d}p}{\mathrm{d}s}(s)\approx \frac{\mathrm{\Delta }p}{\mathrm{\Delta }s}(s)=\frac{p_2(s)-p_1(s)}{l_{\mathrm{t}}}·\end{array}$$(16)

Since the thread has a finite length, gravity changes along this thread; therefore this term is approximated by the average gravity on the thread

$$\begin{array}{cc}\hfill {\overline{g}}_{\Vert }(s)& =\frac{1}{l_{\mathrm{t}}}{{\displaystyle \int }}_{s-l_{\mathrm{t}}/2}^{s+l_{\mathrm{t}}/2}{g}_{\Vert }(s)\mathrm{d}s\hfill \\ \hfill & =-g\frac{2}{\pi }\frac{L}{l_{\mathrm{t}}}sin\left(\frac{\pi }{2}\frac{l_{\mathrm{t}}}{L}\right)cos\left(\pi \frac{s}{L}\right)·\hfill \end{array}$$(17)

Since under typical conditions l_t ≪ L, using the Maclaurin expansion of the sinus for small arguments, we recover exactly g_∥(s) in the previous expression. Hereafter, to simplify things we use g_∥(s) from Eq. (1) instead of Eq. (17).

The hydrostatic equation for the thread is therefore

$$\begin{array}{c}\hfill \frac{p_2(s)-p_1(s)}{l_{\mathrm{t}}}={\rho }_{\mathrm{t}}{g}_{\Vert }(s).\end{array}$$(18)

Using Eq. (3) for the sound speed together with the previous expressions for pressure and gravity we obtain the following equation for s, representing the location under static equilibrium of the dense thread along the tube,

$$\begin{array}{c}\hfill {\left(\frac{L/2-l_{\mathrm{t}}/2}{L-s-l_{\mathrm{t}}/2}\right)}^{\gamma }-{\left(\frac{L/2-l_{\mathrm{t}}/2}{s-l_{\mathrm{t}}/2}\right)}^{\gamma }=-\frac{\gamma }{c_{\mathrm{s}0}^2}\frac{{\rho }_{\mathrm{t}}}{{\rho }_0}g{l}_{\mathrm{t}}cos\left(\pi \frac{s}{L}\right),\end{array}$$(19)

where $c_{\text{s0}}^2=\gamma p_0/{\rho }_0$, and ρ_t/ρ₀ corresponds to the density contrast between the thread and the density background at s = L/2 − l_t/2. The parameter ρ_t/ρ₀ must be much larger than 1 since we are not taking into account any energy flux through the thread and its mass and length are assumed to be constant.

Equation (19) is transcendental and has to be solved numerically. From now on we denote s_e a solution of Eq. (19). It is easy to verify that s_e = L/2 is a solution for any values of the parameters. Hence, there is always an equilibrium solution at the tube apex ($s_{\mathrm{e}}^{\mathcal{A}}=L/2$) that can be stable or unstable. Furthermore, the symmetry of the model with respect to s = L/2 implies that if there are additional solutions, they appear in pairs ($s_{\mathrm{e}}^{\mathcal{L}},s_{\mathrm{e}}^{\mathcal{R}}$) verifying $s_{\mathrm{e}}^{\mathcal{L}}+s_{\mathrm{e}}^{\mathcal{R}}=L$, to wit, the solutions are distributed symmetrically with respect to s = L/2: one on the left ($s_{\mathrm{e}}^{\mathcal{L}}$) and the other on the right side ($s_{\mathrm{e}}^{\mathcal{R}}$) of the tube.

To understand the nature of these possible situations, in Fig. 2 we plotted the left- and right-hand side terms of Eq. (19) as a function of s, for fixed parameters L, c_s0, g, ρ_t/ρ₀, and varying l_t using the projected gravity given in Eq. (1). The gravity acceleration is g = 274 m s⁻², and we take a spatial reference length of H = 10 000 km. The sound speed takes a value of c_s0 = 166 km s⁻¹ for a fully ionised coronal tube with 10⁶ K temperature. When the two curves cross (see circles) we have a solution at that particular point. On the top panel we have only one solution: the trivial one at $s_{\mathrm{e}}^{\mathcal{A}}=L/2$. On the bottom panel however three solutions appear: the trivial one plus two complementary solutions. As we describe in Sect. 4, in the first situation the mass remains at its original position, and, upon being subject to a small disturbance, it oscillates around this point because this equilibrium is stable. In the second situation, the mass is not stable at the initial position at the apex, and after being perturbed the mass moves to one of the solutions along the legs that are energetically more favourable. The stability analysis performed in this and in the following sections provides a simple criterion that tells us when a plasma thread is gravitationally stable or unstable.

The transition between one to three solutions is represented in the middle panel of Fig. 2 and the mathematical condition is that the slopes of the curves are the same for the left- and right-hand sides of Eq. (19),

Fig. 2. Left- (black line) and right- (red line) hand sides of Eq. (19) as a function of s/L, for three different values of lt. Top panel: system is linearly stable with one solution (see circle). Bottom panel: system relaxes to a new equilibrium at a lower height (see the two complementary solutions) and the trivial solution becomes an unstable solution. In all cases, L = 10 H and ρt = 100 ρ0; lt = L/110 (top panel), lt = ltb (middle panel), and lt = L/30 (bottom panel). The value ltb ≈ L/77.08 is obtained when the derivatives of the two curves at s = L/2 are the same (see Eq. (20)). The vertical dashed lines represent the edges of the domain of the variable s, namely s = lt/2 and s = L − lt/2.

$$\begin{array}{c}\hfill \frac{{(L/2-l_{\mathrm{t}}/2)}^{\gamma }}{{(L-s-l_{\mathrm{t}}/2)}^{\gamma +1}}+\frac{{(L/2-l_{\mathrm{t}}/2)}^{\gamma }}{{(s-l_{\mathrm{t}}/2)}^{\gamma +1}}=\frac{1}{c_{\mathrm{s}0}^2}\frac{{\rho }_{\mathrm{t}}}{{\rho }_0}g{l}_{\mathrm{t}}\frac{\pi }{L}sin\left(\pi \frac{s}{L}\right)·\end{array}$$(20)

Considering the following equation:

$$\begin{array}{c}\hfill f(s)=\frac{\mathrm{d}}{\mathrm{d}s}[-\frac{1}{{\rho }_{\mathrm{t}}}\frac{\mathrm{d}p}{\mathrm{d}s}(s)+g_{\Vert }(s)]=0,\end{array}$$(21)

and using the expressions for the pressure derivative (Eq. (16)) and the projected gravity (Eq. (1)), we obtain Eq. (20). When Eq. (20) or Eq. (21) is satisfied together with Eq. (19) the system is said to have a bifurcation point. We realise that the left-hand term of Eq. (21) is precisely the term inside the square root of the computed eigenvalues in Eqs. (12) and (13). We have shown that depending on the sign of f(s_e) the nature of the equilibrium solutions changes. When f(s_e) is positive (Eq. (12)) the eigenvalue λ₋ represents an exponentially damped solution that is eventually dominated by the growing solution in time associated to λ₊ of the form e^(τt) where the growth time is τ = λ₊. When f(s_e) is negative we have a stable linear centre with a harmonic time-dependence of the form e^(iωt) and frequency ω = λ₊/i (Eq. (13)). In this case, it is known from linear theory of ODE systems that we cannot conclude that the equilibrium point is truly a centre since it could be also a stable or unstable spiral (e.g., Jordan & Smith 1987). Nevertheless, when the system is reversible or energy conservative (Hamiltonian), and in our case the two conditions are fulfilled, the linear centre is a true centre (e.g., Strogatz 2018 Sects. 6.5 and 6.6) and stability is ensured. In Sect. 4 we confirm, using physical arguments, that the equilibrium point is stable.

We discuss the results from the physical point of view. Figure 2 shows that the existence of the non-trivial solutions of Eq. (19) depends on the thread model parameters. Figure 3 (top panel) shows the equilibrium position s_e as a function of ρ_t/ρ₀ for different values of L when the thread length is l_t = 0.1 H. We obtain that the system only develops the $s_{\mathrm{e}}^{\mathcal{L}}$ and the $s_{\mathrm{e}}^{\mathcal{R}}$ solutions for values of ρ_t beyond ρ_tb. When ρ_t < ρ_tb, meaning before the appearance of the bifurcation, $s_{\mathrm{e}}^{\mathcal{A}}$ is the only equilibrium position, so that we are in the situation of the top panel of Fig. 2 and therefore $s_{\mathrm{e}}^{\mathcal{A}}/L=0.5$ is stable. However, when the thread structure enables the existence of two complementary solutions, these are energetically more suitable so that $s_{\mathrm{e}}^{\mathcal{A}}$, which remains located at s_e/L = 0.5 (see horizontal dashed lines in Fig. 3 top panel), becomes an unstable solution. To distinguish between stable or unstable solutions, we represent in Fig. 3 (and in the following figures) stable solutions as solid curves and unstable solutions as dashed curves. In addition, we realise that for the three values of L, $s_{\mathrm{e}}^{\mathcal{L}}$ and $s_{\mathrm{e}}^{\mathcal{R}}$ appear around the same bifurcation value ρ_tb. This means that the bifurcation density value from which we obtain an equilibrium out of the apex of the structure, does not strongly depend on the length of the magnetic tube. We note that s_e is normalised to H for normalisation purposes since if we keep the normalisation to L, the solutions of the three configurations overlap. This means that the normalised s_e/L does not depend strongly on L for this set of parameters. This feature is clearer in the bottom panel of Fig. 3, where s_e as a function of L is plotted; when L is large the stable solutions are almost horizontal. Moreover, as was expected, the equilibrium position of the thread drops with increasing ρ_t. It is interesting to study the dependence of the equilibrium position on the thread length. Figure 3 (middle panel) shows s_e as a function of l_t/L for different values of L when ρ_t/ρ₀ = 100. For L = 10 H we obtain that below the bifurcation $s_{\mathrm{e}}=s_{\mathrm{e}}^{\mathcal{A}}$ is the only solution of Eq. (19) and this is stable, but when l_t > l_tb, $s_{\mathrm{e}}^{\mathcal{A}}$ becomes unstable and the emerged branches $s_{\mathrm{e}}^{\mathcal{L}}$ and $s_{\mathrm{e}}^{\mathcal{R}}$ are stable. For L = 15 H and L = 20 H the bifurcation thread length is out of the plot and their $s_{\mathrm{e}}^{\mathcal{A}}$ overlaps with that for L = 10 H at s_e/L = 0.5. In Fig. 3 (middle panel) we also see that the equilibrium position of the thread changes along the magnetic tube when the thread length increases. Figure 3 (top and middle panels) show the same behaviour for all sets of parameters, which consists in a sole stable solution for values below the bifurcation, and two symmetric stable solutions out of the apex plus an unstable solution at s_e = L/2 beyond the bifurcation point. However, Fig. 3 (bottom panel) shows that for some set of parameters, the system does not present the $s_{\mathrm{e}}^{\mathcal{L}}$ and $s_{\mathrm{e}}^{\mathcal{R}}$ equilibria when L varies so that $s_{\mathrm{e}}^{\mathcal{A}}$ is always stable. Interestingly, the bottom panel of Fig. 3 shows that for small changes of ρ_t, $s_{\mathrm{e}}^{\mathcal{L}}$ and $s_{\mathrm{e}}^{\mathcal{R}}$ vary significantly. Moreover, we see that the bifurcation length L_b for ρ_t/ρ₀ = 133 and ρ_t/ρ₀ = 138 is short, around L_b = 2.7 H and L_b = 1.4 H, respectively. These values are shorter than the typical magnetic tube observations.

Fig. 3. Solutions se of Eq. (19) as a function of ρt/ρ0 (top panel), lt/L (middle panel), and L/H (bottom panel). The expression lt = 0.1 H in the top and bottom panels and ρt/ρ0 = 100 in the middle panel. The solid lines correspond to stable solutions and dashed lines to unstable solutions. The thin vertical dashed lines correspond to the bifurcation values that have been calculated by solving Eqs. (19) and (20) simultaneously and coincide with the results in Eq. (32) (top panel), Eq. (33) (middle panel), and Eq. (34) (bottom panel).

More aspects about the bifurcation values and how to calculate them are given in Sect. 4. But we finish this section by analysing the presented results from a dynamical systems perspective. The transition from one stable equilibrium solution to two symmetric stable solutions and one unstable solution is characteristic of the supercritical pitchfork bifurcation (see e.g., Sect. 3.4 of Strogatz 2018 or Sect. 20.1E of Wiggins 2003). This bifurcation is present in systems that have a symmetry, such as the symmetry between the left and right sides of our magnetic tube when the thread is initially at the tube apex. In Fig. 3 we find the same type of pitchfork bifurcation in the three panels, although in the bottom panel the unstable solutions for ρ_t/ρ₀ = 133 and ρ_t/ρ₀ = 138 are overplotted by the stable solution of ρ_t/ρ₀ = 129.

3.2. Thread initially at any location along the tube

It is straightforward to derive analytical expressions for the situation of a thread initially located at any position along the tube, denoted by s₀, and not necessarily at L/2 as discussed in the previous section. We simply need to rewrite Eq. (15) as

$$\begin{array}{cc}& p_1(s)=p_0{\left(\frac{s_0-l_{\mathrm{t}}/2}{s-l_{\mathrm{t}}/2}\right)}^{\gamma },\hfill \\ \hfill & p_2(s)=p_0{\left(\frac{L-s_0-l_{\mathrm{t}}/2}{L-s-l_{\mathrm{t}}/2}\right)}^{\gamma },\hfill \end{array}$$(22)

and the corresponding equilibrium equation analogous to Eq. (19) becomes

$$\begin{array}{c}\hfill {\left(\frac{L-s_0-l_{\mathrm{t}}/2}{L-s-l_{\mathrm{t}}/2}\right)}^{\gamma }-{\left(\frac{s_0-l_{\mathrm{t}}/2}{s-l_{\mathrm{t}}/2}\right)}^{\gamma }=-\frac{\gamma }{c_{\mathrm{s}0}^2}\frac{{\rho }_{\mathrm{t}}}{{\rho }_0}g{l}_{\mathrm{t}}cos\left(\pi \frac{s}{L}\right)·\end{array}$$(23)

Equation (23) simplifies to Eq. (19) when s₀ = L/2. Different types of representative solutions are shown in Fig. 4 and in this case the parameter that has been chosen to change is ρ_t (in Fig. 2 it was l_t). In these plots the initial position of the thread is s₀ = 0.4/L, therefore located on the left half of the tube. The left-hand side of Eq. (23) is not symmetric anymore with respect to the tube apex, and this facilitates the existence either of one, two, or three solutions, according to Fig. 4. But we note that by continuity there must be a situation in these circumstances with two solutions when the left- and right-hand curves are tangential. This case is calculated as in Sect. 3.1 by imposing that the spatial derivatives of the left- and right-hand sides of Eq. (23) are equal as follows:

$$\begin{array}{c}\hfill \frac{{(L-s_0-l_{\mathrm{t}}/2)}^{\gamma }}{{(L-s-l_{\mathrm{t}}/2)}^{\gamma +1}}+\frac{{(s_0-l_{\mathrm{t}}/2)}^{\gamma }}{{(s-l_{\mathrm{t}}/2)}^{\gamma +1}}=\frac{1}{c_{\mathrm{s}0}^2}\frac{{\rho }_{\mathrm{t}}}{{\rho }_0}g{l}_{\mathrm{t}}\frac{\pi }{L}sin\left(\pi \frac{s}{L}\right)·\end{array}$$(24)

Fig. 4. Left- (black line) and right- (red line) hand sides of Eq. (23) as a function of s/L, for three different situations with one, two, and three solutions. In these plots s0 = 0.4/L (where the black line crosses zero), L = 10 H, lt = L/50; ρt/ρ0 = 100 (top panel), ρt/ρ0 = ρtb/ρ0 (middle panel), and ρt/ρ0 = 200 (bottom panel). The value ρtb/ρ0 ≈ 137.28 is obtained from Eqs. (23) and (24). The vertical dashed lines represent the edges of the domain of the variable s, namely s = lt/2 and s = L − lt/2.

Thus, Eq. (24) has to be satisfied together with Eq. (23) for the system to have two equilibrium solutions. An example of such a situation is found in Fig. 4 middle panel. For simplicity, we still denote the equilibrium position on the left side of the tube as $s_{\mathrm{e}}^{\mathcal{L}}$, as $s_{\mathrm{e}}^{\mathcal{R}}$ on the right side, and as $s_{\mathrm{e}}^{\mathcal{A}}$ for the central near apex solution, so that when the system has only two equilibrium solutions, $s_{\mathrm{e}}^{\mathcal{R}}=s_{\mathrm{e}}^{\mathcal{A}}$ (or $s_{\mathrm{e}}^{\mathcal{L}}=s_{\mathrm{e}}^{\mathcal{A}}$). As in Sect. 3.1 the sign of the difference between the left-hand term and the right-hand term of Eq. (24) reveals the nature of the different solutions.

In the bifurcation diagram we find that the characteristic pitchfork for the symmetric case disconnects into two curves: one below s₀/L = 0.5 and another above this value. An example is shown in Fig. 5 (top panel) when the density contrast is changed. The lower curve is always stable, whereas the upper curve has both stable and unstable branches. Therefore, the pitchfork bifurcation occurring in the symmetrical case is essentially replaced in the non-symmetrical case by a stable state (lower curve) alongside a saddle-node bifurcation (upper curve). As we increase the density contrast there is no longer a sharp transition at the bifurcation point representative of the symmetric case (compare with Fig. 3). Hence, the equilibrium position of the thread varies smoothly with the density contrast on the lower stable curve. Furthermore, the stable branch of the upper curve, which is also a permitted equilibrium state of the system, is not accessible unless the thread is given a large initial disturbance, thereby making this equilibrium more favourable. Qualitatively similar results are obtained for the bifurcation diagrams as a function of the thread length and the total length in the middle and bottom panels of Fig. 5, respectively. As was expected, when s₀/L → 0.5, the diagrams approach those in Fig. 3, characterised by a stable solution at s_e/L = 0.5 below the bifurcation values; however we see that the two curves (upper and lower) still remain disconnected for s₀/L = 0.49.

Fig. 5. Same as Fig. 3 but for Eq. (23) and for different values of s0/L. L = 10 H and lt = 0.2 H (top panel); ρt/ρ0 = 100 and L = 10 H (middle panel); lt = 0.2 H and ρt/ρ0 = 137 (bottom panel). The bifurcation points (intersections between vertical and horizontal thin dashed lines) have been calculated by solving Eqs. (23) and (24) simultaneously.

Figure 6 shows the situation when there is a transition between one to three solutions (orange lines, ρ_t/ρ₀ = 200) when the initial position of the thread, s₀, is changed. One of the emerging solutions around s₀/L = 0.3, $s_{\mathrm{e}}^{\mathcal{A}}$, is unstable, while the other solution, $s_{\mathrm{e}}^{\mathcal{R}}$, is stable. This behaviour is similar to that found for the symmetric case with s₀ = L/2 described earlier. In Fig. 6 a case with no bifurcations but with multiple solutions for fixed parameter values (ρ_t/ρ₀ = 400) is also represented. The solutions that are around the apex (dashed curves) have a different behaviour with respect to the left- and right-hand solutions. In this case $s_{\mathrm{e}}^{\mathcal{A}}$ decreases when s₀ is raised, while for the other two solutions $s_{\mathrm{e}}^{\mathcal{L}}$ and $s_{\mathrm{e}}^{\mathcal{R}}$ increase with s₀. When s₀ moves closer to the tube apex, $s_{\mathrm{e}}^{\mathcal{A}}\to L/2$ and this solution is unstable. In conclusion, given that the initial thread position, s₀, is not at the apex, this parameter breaks the symmetry of the system. In bifurcation theory it is called an imperfection parameter. The resulting bifurcation diagram of Fig. 6 is no longer symmetric about a horizontal line and possesses an S-shape (see Wang 1994 and Sect. 3.6 of Strogatz 2018). The orange curve in this figure has a lower and an upper bifurcation point, s_0, l/L ≃ 0.3 and s_0, u/L ≃ 0.7, respectively. This S-shaped bifurcation curve can be thought of as the composition of two saddle-node bifurcation curves, the lower (upper) one with its bifurcation point at s₀ = s_0, l (s₀ = s_0, u). The black curve of Fig. 6 is a special case of an S-shaped bifurcation curve because the lower and upper bifurcation points are outside the range of validity of the parameter s₀. Therefore, the system always supports three solutions.

Fig. 6. Solutions se of Eq. (23) as a function of s0/L for ρt/ρ0 = 200 (orange lines) and ρt/ρ0 = 400 (black lines). The expression L = 10 H and lt = 0.2 H. The curves correspond to $s_{\mathrm{e}}^{\mathcal{L}}$ and $s_{\mathrm{e}}^{\mathcal{R}}$ (solid lines) and $s_{\mathrm{e}}^{\mathcal{A}}$ (dashed line). The bifurcation points around s0/L = 0.3 and s0/L = 0.7 (intersections between vertical and horizontal orange dashed lines) have been calculated by solving Eqs. (23) and (24) simultaneously.

4. Linear stability analysis

We have intuitively described the nature of stable and unstable solutions in the previous section and complemented this description with known results from bifurcation theory in mathematics. In the present section we provide a precise stability analysis based on linear theory but using physical grounds. We derive an expression for the oscillatory frequency of the thread, ω, as a function of the equilibrium parameters. The sign of ω² renders estimable knowledge about the stability properties of the configuration.

4.1. Thread initially at the tube apex

We start with the trivial solution at the tube apex. We consider a small longitudinal displacement, δs, of the centre of the plasma thread. We assume that δs ≪ L because we are interested in the linear regime. The restoring forces acting on the plasma thread are the pressure gradient and the projected gravity force. Let us write the pressure at both sides of the thread as a function of the displacement with respect to L/2, the equilibrium position. Using the adiabatic assumption again, the pressure on the left and right parts of the tube is given by

$$\begin{array}{cc}& p_1(\delta s)=p_0{\left(\frac{L/2-l_{\mathrm{t}}/2}{L/2-l_{\mathrm{t}}/2+\delta s}\right)}^{\gamma },\hfill \\ \hfill & p_2(\delta s)=p_0{\left(\frac{L/2-l_{\mathrm{t}}/2}{L/2-l_{\mathrm{t}}/2-\delta s}\right)}^{\gamma }.\hfill \end{array}$$(25)

The terms in parenthesis are written as

$$\begin{array}{cc}\hfill {\left(\frac{L/2-l_{\mathrm{t}}/2}{L/2-l_{\mathrm{t}}/2\pm \delta s}\right)}^{\gamma }& ={\left(\frac{1}{1\pm 2\delta s/(L-l_{\mathrm{t}})}\right)}^{\gamma }\hfill \\ \hfill & ={\left(1\pm 2\delta s/(L-l_{\mathrm{t}})\right)}^{-\gamma }\hfill \\ \hfill & \simeq 1\mp 2\gamma \delta s/(L-l_{\mathrm{t}}),\hfill \end{array}$$(26)

where in the last step we used the Maclaurin series approximation for small arguments (because δs/L ≪ 1). This approximation is crucial to find the linear result we are looking for, since pressure fluctuations are now linearly proportional to δs. We have that

$$\begin{array}{cc}& p_1(\delta s)\simeq p_0[1-2\gamma \delta s/(L-l_{\mathrm{t}})],\hfill \\ \hfill & p_2(\delta s)\simeq p_0[1+2\gamma \delta s/(L-l_{\mathrm{t}})],\hfill \end{array}$$(27)

and the pressure gradient is given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}p}{\mathrm{d}s}(\delta s)\simeq \frac{\mathrm{\Delta }p}{\mathrm{\Delta }s}(\delta s)=\frac{p_2-p_1}{l_{\mathrm{t}}}=\frac{4p_0\gamma \delta s}{l_{\mathrm{t}}(L-l_{\mathrm{t}})}·\end{array}$$(28)

For the projected gravity around s = L/2 it is not difficult to see that Eq. (1) reduces to a sinus of δs,

$$\begin{array}{c}\hfill g_{\Vert }(L/2+\delta s)=gsin\left(\pi \frac{\delta s}{L}\right)\simeq g\frac{\pi }{L}\delta s,\end{array}$$(29)

where we again used the approximation for small arguments in the last step.

We have expressions for all the terms that appear in the momentum equation, that is Eq. (6) for s = L/2 + δs, which is now

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{t}}\frac{{\mathrm{d}}^2\delta s}{\mathrm{d}{t}^2}& =-\frac{4p_0\gamma }{l_{\mathrm{t}}(L-l_{\mathrm{t}})}\delta s+{\rho }_{\mathrm{t}}g\frac{\pi }{L}\delta s\hfill \\ \hfill & =\delta s[-\frac{4p_0\gamma }{l_{\mathrm{t}}(L-l_{\mathrm{t}})}+{\rho }_{\mathrm{t}}g\frac{\pi }{L}]·\hfill \end{array}$$(30)

The right-hand side term of this equation is proportional to δs only. The solution to the differential equation is of the form e^iωt and we have that

$$\begin{array}{c}\hfill {\omega }^2=\frac{4c_{\mathrm{s}0}^2}{l_{\mathrm{t}}(L-l_{\mathrm{t}}){\rho }_{\mathrm{t}}/{\rho }_0}-g\frac{\pi }{L},\end{array}$$(31)

where the sound speed of the background was used. According to Eq. (31) we realise that we can have the following situations: ω² > 0 (ω real), meaning that the motion is purely oscillatory around a stable equilibrium; ω² < 0 (ω purely imaginary), indicating that the system is unstable; and ω² = 0. This is the expected situation based on the results of Sect. 3 about the equilibrium. Interestingly, it is not difficult to see that the growth rates and frequencies derived from the dispersion relation in Eq. (31) are exactly the same τ and ω found in Sect. 3.1 and derived from a mathematical point of view.

The transition between the two regimes (ω² = 0) provides helpful information because it corresponds to a bifurcation point, where, as we have seen, the number of solutions and/or their stability can change. This transition takes place when the following condition for the density contrast is satisfied:

$$\begin{array}{c}\hfill {{\rho }_{\mathrm{t}}}_{\mathrm{b}}/{\rho }_0=\frac{4c_{\mathrm{s}0}^{2}L}{g\pi l_{\mathrm{t}}(L-l_{\mathrm{t}})},\end{array}$$(32)

where the density contrasts below this bifurcation value correspond to a stable solution, whereas for ρ_t > ρ_tb the stable solution turns unstable. This behaviour is shown in Fig. 7 (top panel), in which ω² is plotted as a function of ρ_t/ρ₀ for different values of L and for the three possible equilibrium positions $s_{\mathrm{e}}^{\mathcal{A}}$, $s_{\mathrm{e}}^{\mathcal{L}}$, and $s_{\mathrm{e}}^{\mathcal{R}}$. Since the initial position of the thread is at the apex of the structure, $s_{\mathrm{e}}^{\mathcal{L}}$ and $s_{\mathrm{e}}^{\mathcal{R}}$ are symmetric and have the same frequency, so only one solid curve with ω² > 0 is visible to the right of the bifurcation point. Frequencies associated with $s_{\mathrm{e}}^{\mathcal{A}}$ are found at the right of the bifurcation point and below zero in Fig. 7. We see that the bifurcation density, which delimits stable solutions (solid curves) with unstable (dashed curves), does not strongly depend on L. This behaviour agrees with Eq. (32), which simplifies to ${{\rho }_{\mathrm{t}}}_{\mathrm{b}}/{\rho }_0=\frac{4c_{\mathrm{s}0}^2}{g\pi l_{\mathrm{t}}}$ when l_t is neglected with respect to L.

Fig. 7. ω2 as a function of ρt/ρ0 (top panel), lt/L (middle panel), and L/H (bottom panel) for $s_{\mathrm{e}}=s_{\mathrm{e}}^{\mathcal{A}}$ (below bifurcation points and lower branch curves, Eq. (31)) and $s_{\mathrm{e}}=s_{\mathrm{e}}^{\mathcal{L}}=s_{\mathrm{e}}^{\mathcal{R}}$ (upper branch curves; Eq. (40)). In these plots s0 = L/2. lt = 0.1 H in the top and bottom panels and ρt/ρ0 = 100 in the middle panel. The bifurcation values (thin dashed vertical lines) have been calculated using the analytical expressions given by Eqs. (32)–(34).

Assuming that ρ_t/ρ₀ and L are fixed, then we have a bifurcation length for the thread that determines if the solution is stable or unstable as follows:

$$\begin{array}{c}\hfill {l_{\mathrm{t}}}_{\mathrm{b}}=\frac{L}{2}(1-\sqrt{1-\frac{16c_{\mathrm{s}0}^2}{g\pi L{\rho }_{\mathrm{t}}/{\rho }_0}}),\end{array}$$(33)

(see middle panel of Fig. 7). For typical prominence values the term inside the square root is always positive and the square root of this term is smaller than 1, providing a physically acceptable value (l_tb > 0). The equilibrium positions of threads with lengths shorter than l_tb are always stable.

Finally, from Eq. (31) we find that the bifurcation length of the magnetic field lines is given by

$$\begin{array}{c}\hfill L_{\mathrm{b}}=\frac{l_{\mathrm{t}}}{1-\frac{4c_{\mathrm{s}0}^2}{g\pi l_{\mathrm{t}}{\rho }_{\mathrm{t}}/{\rho }_0}},\end{array}$$(34)

and the term in the denominator is not necessarily positive. A negative bifurcation length is non-physical, meaning that in this situation it is not possible to find a transition from the stable to the unstable regime. Therefore, when

$$\begin{array}{c}\hfill \frac{4c_{\mathrm{s}0}^2}{g\pi l_{\mathrm{t}}{\rho }_{\mathrm{t}}/{\rho }_0}>1,\end{array}$$(35)

the bifurcation length is negative, thereby implying in our case that the solution is always stable (ω² > 0) independent of the length L; this is easy to check introducing the condition given by Eq. (35) into Eq. (31). Figure 7 (bottom panel) shows ω² as a function of L for different values of ρ_t and we find a case in which Eq. (35) is satisfied and ω² is always positive. For values of the parameters that make the left-hand side of Eq. (35) less than 1, the solution is either stable or unstable, depending on the length of the field lines, if L < L_b the system that only develops $s_{\mathrm{e}}^{\mathcal{A}}$ is stable and if L > L_b it is unstable (see lower branch curves in bottom panel of Fig. 7). Interestingly, the bifurcation length L_b depends strongly on the density contrast (see the short range of ρ_t/ρ₀ in Fig. 7 bottom panel). We note that to have a stable solution the three criteria must be simultaneously satisfied; these criteria are ρ_t < ρ_tb, l_t < l_tb, and L < L_b or equivalently Eq. (35) for the length of the field lines. This is an interesting result that provides useful information about the behaviour of the system depending on the parameters.

It is worthy noting that the first term on the right-hand side of Eq. (31) is the frequency squared of a slow mode in a slab model of length l_t inside a field line with a length of L (Joarder & Roberts 1992; Díaz et al. 2010; Soler et al. 2010). Luna et al. (2012) derived a similar expression for the frequency but for a magnetic field with a dip at the location of the thread. In their formula the sign in front of the gravitational term is positive and therefore the instability is not present in their configuration; their system is always (linearly) stable because they consider a concave upwards geometry for the magnetic field.

Now we generalise the frequency of oscillation for any equilibrium position, s_e, along the tube and hence we do not restrict the analysis to the equilibrium solution at the tube apex. But we note that the thread is still initially located at s = L/2. Under these conditions we have that

$$\begin{array}{cc}& p_1(\delta s)=p_{10}{\left(\frac{s_{\mathrm{e}}-l_{\mathrm{t}}/2}{s_{\mathrm{e}}-l_{\mathrm{t}}/2+\delta s}\right)}^{\gamma },\hfill \\ \hfill & p_2(\delta s)=p_{20}{\left(\frac{s_{\mathrm{e}}-l_{\mathrm{t}}/2}{L-s_{\mathrm{e}}-l_{\mathrm{t}}/2-\delta s}\right)}^{\gamma },\hfill \end{array}$$(36)

where p₁₀ and p₂₀ correspond to the equilibrium pressures that satisfy Eq. (15). We rewrite p₁ as

$$\begin{array}{c}\hfill p_1(\delta s)\simeq p_{10}[1-\gamma \delta s/(s_{\mathrm{e}}-l_{\mathrm{t}}/2)],\end{array}$$(37)

where we used the Maclaurin series approximation for small arguments, but if s_e is similar to l_t/2 the approximation is not valid. An analogous expression is derived for p₂. The pressure gradient using the previous equations is approximated by

$$\begin{array}{cc}\hfill \frac{\mathrm{d}p}{\mathrm{d}s}(\delta s)\simeq & \frac{p_2(\delta s)-p_1(\delta s)}{l_{\mathrm{t}}}\hfill \\ \hfill & =\frac{1}{l_{\mathrm{t}}}(p_{20}-p_{10}+\frac{p_{20}\gamma \delta s}{L-s_{\mathrm{e}}-l_{\mathrm{t}}/2}+\frac{p_{10}\gamma \delta s}{s_{\mathrm{e}}-l_{\mathrm{t}}/2})·\hfill \end{array}$$(38)

For the projected gravity we need to perform a Taylor expansion of the cosine around the equilibrium position, which is now s_e,

$$\begin{array}{c}\hfill cos\left(\pi \frac{s}{L}\right)=cos\left(\pi \frac{s_{\mathrm{e}}}{L}\right)-\frac{\pi }{L}sin\left(\pi \frac{s_{\mathrm{e}}}{L}\right)(s-s_{\mathrm{e}})+\dots ,\end{array}$$(39)

where we only retain up to linear terms in s − s_e. We have that δs = s − s_e according to our notation.

In the momentum equation, using Eqs. (38) and (39), we find two terms that are not proportional to δs but that cancel out because they correspond to the equilibrium condition given by Eq. (19). We derive that the square of the frequency is given by

$$\begin{array}{cc}\hfill {\omega }^2=& \frac{c_{\mathrm{s}0}^2}{l_{\mathrm{t}}{\rho }_{\mathrm{t}}/{\rho }_0}[\frac{{(L/2-l_{\mathrm{t}}/2)}^{\gamma }}{{(L-s_{\mathrm{e}}-l_{\mathrm{t}}/2)}^{\gamma +1}}+\frac{{(L/2-l_{\mathrm{t}}/2)}^{\gamma }}{{(s_{\mathrm{e}}-l_{\mathrm{t}}/2)}^{\gamma +1}}]\hfill \\ \hfill & -g\frac{\pi }{L}sin\left(\pi \frac{s_{\mathrm{e}}}{L}\right)·\hfill \end{array}$$(40)

When s_e = L/2, that is the trivial solution, Eq. (40) is exactly the same as Eq. (31), recovering the previous situation. The solutions of Eq. (40) are shown in Fig. 7. Below the bifurcation values ρ_tb, l_tb, and L_b, the equilibrium solution is $s_{\mathrm{e}}^{\mathcal{A}}$. When the bifurcation values are reached, $s_{\mathrm{e}}^{\mathcal{A}}$ becomes unstable (ω² < 0), and as expected, $s_{\mathrm{e}}^{\mathcal{L}}$ is stable (ω² > 0). For $s_{\mathrm{e}}^{\mathcal{R}}$ the value of frequency is exactly the same as that for $s_{\mathrm{e}}^{\mathcal{L}}$ and therefore only one curve is visible in the plot. For $s_{\mathrm{e}}^{\mathcal{L}}$ and $s_{\mathrm{e}}^{\mathcal{R}}$ (upper branch curves) the frequency squared monotonically increases with ρ_t and l_t, but it increases with L until a maximum and then it starts decreasing with L. It is important to mention that all the results shown in Fig. 7 are the solutions of the configurations shown in Fig. 3 and the stability analysis done in Sect. 3.1 agrees with the results in the present section.

4.2. Thread initially at any location along the tube

Finally, using the same procedure as before we derive an expression for the oscillation frequency of a thread around the final equilibrium position s_e, which has been achieved from an initial location of the thread at s₀ ≠ L/2. The corresponding pressures at the edge of the thread due to a small displacement δs from the equilibrium position are expressed as Eq. (36), but now p₁₀ and p₂₀ correspond to the equilibrium pressures that satisfy Eq. (22). Using similar approximations to those in Eq. (37) and cancelling the equilibrium terms according to Eq. (23), we eventually find that the frequency squared is written as

$$\begin{array}{cc}\hfill {\omega }^2=& \frac{c_{\mathrm{s}0}^2}{l_{\mathrm{t}}{\rho }_{\mathrm{t}}/{\rho }_0}[\frac{{(L-s_0-l_{\mathrm{t}}/2)}^{\gamma }}{{(L-s_{\mathrm{e}}-l_{\mathrm{t}}/2)}^{\gamma +1}}+\frac{{(s_0-l_{\mathrm{t}}/2)}^{\gamma }}{{(s_{\mathrm{e}}-l_{\mathrm{t}}/2)}^{\gamma +1}}]\hfill \\ \hfill & -g\frac{\pi }{L}sin\left(\pi \frac{s_{\mathrm{e}}}{L}\right)·\hfill \end{array}$$(41)

The previous frequency expressions given in Eq. (31) when s₀ = s_e = L/2 and Eq. (40) when s₀ = L/2, s_e ≠ L/2 are recovered in Eq. (41). In this equation s_e must be a solution to Eq. (23). Obtaining analytical bifurcation values, except for the trivial case (s₀ = s_e = L/2), is much more difficult now and it is only possible for the bifurcation density contrast. The rest of bifurcation values are determined numerically. As expected from bifurcation theory, imposing the condition for a bifurcation value, that is ω² = 0 in Eq. (41), is completely equivalent to the condition given by Eq. (24) that determines the bifurcation point based only on the equilibrium configuration. This is confirmed by the comparison of the bifurcation values in Fig. 8 and the respective results in Fig. 5. This is also a solid indication that the performed eigenfrequency analysis in this section is coherent.

Fig. 8. ω2 as function of different equilibrium parameters (Eq. (41)) for three different initial positions of the thread, s0. The expression L = 10 H and lt = 0.2 H (top panel), ρt/ρ0 = 100 and L = 10 H (middle panel), and lt = 0.2 H and ρt/ρ0 = 137 (bottom panel). The vertical lines correspond to the bifurcation points and are calculated by imposing that ω2 = 0.

Now the break of symmetry due to the imperfection parameter s₀ causes the difference between the frequency squared of $s_{\mathrm{e}}^{\mathcal{L}}$ and $s_{\mathrm{e}}^{\mathcal{R}}$. This is visible in the plotted upper curves of Fig. 8 (solid curves). Beyond the bifurcation values two branches appear that correspond to $s_{\mathrm{e}}^{\mathcal{R}}$ and $s_{\mathrm{e}}^{\mathcal{A}}$, which are positive for $s_{\mathrm{e}}^{\mathcal{R}}$ (solid curves) and negative for $s_{\mathrm{e}}^{\mathcal{A}}$ (dashed curves). Logically, when s₀/L → 0.5 (purple curves), ω² for $s_{\mathrm{e}}^{\mathcal{R}}$ approaches the same value as for $s_{\mathrm{e}}^{\mathcal{L}}$. For $s_{\mathrm{e}}^{\mathcal{R}}$ (emerging at the bifurcation point), ω² increases with ρ_t/ρ₀ (top panel) and l_t/H (middle panel). For $s_{\mathrm{e}}^{\mathcal{L}}$, stable in the whole parameter range, ω² decreases with ρ_t/ρ₀ and l_t/H until it reaches a minimum located before the bifurcation point from where it starts increasing. On the contrary, in the bottom panel of Fig. 8 we see that for $s_{\mathrm{e}}^{\mathcal{L}}$ the frequency decreases for large values of L.

In Fig. 9 the three square of the frequencies corresponding to the solutions in Fig. 6 are plotted as a function of the initial position of the thread, s₀. To calculate these frequencies, we used the values of s_e that satisfy Eq. (23). We note that $s_{\mathrm{e}}^{\mathcal{L}}$ and $s_{\mathrm{e}}^{\mathcal{R}}$ are always stable while $s_{\mathrm{e}}^{\mathcal{A}}$ is always unstable. As expected from physical grounds the left and right side solutions converge to the same value when there is symmetry in the system, that is for s₀ = L/2; we note that beyond this point the meaning of $s_{\mathrm{e}}^{\mathcal{L}}$ and $s_{\mathrm{e}}^{\mathcal{R}}$ has to be inverted. We also observe the transition between one and three solutions for a specific choice of parameters (ρ_t/ρ₀ = 200, orange curves). At the bifurcation point two solutions emerge from ω = 0: one representing a stable solution ($s_{\mathrm{e}}^{\mathcal{R}}$) and another corresponding to an unstable solution ($s_{\mathrm{e}}^{\mathcal{A}}$). We also see that ω² significantly varies with s₀/L for stable solutions, but it is essentially independent of s₀/L for unstable solutions.

Fig. 9. ω2 as a function of s0/L (Eq. (41)) for the three solutions obtained in Fig. 6. The solid lines correspond to $s_{\mathrm{e}}^{\mathcal{L}}$ (lines decreasing from left to right) and $s_{\mathrm{e}}^{\mathcal{R}}$ (lines increasing from left to right), and the dashed lines to $s_{\mathrm{e}}^{\mathcal{A}}$. The horizontal line represents ω = 0 and divides the plot into stable and unstable solutions. In this plot the left and right solutions are stable while the near apex solution is unstable. The bifurcation points (vertical dashed lines) have been calculated by solving Eqs. (23) and (24) simultaneously, and it is equivalent to impose ω2 = 0 in Eq. (41) once se is known.

5. Conclusions and discussion

The thread-like structures studied in this work appear to be present in filaments everywhere (Lin et al. 2008), along the horizontal spines or the quasi-vertical barbs. High-resolution observations of the solar corona show that the fine structure of prominences is very dynamic. Even during the prominence formation we observe plasma flows along the magnetic fields. The localised heating above the flux tube footpoints produces a cyclic pattern of evaporation of the chromospheric plasma, which condenses in the coronal part of the tube and produces the plasma flows. During the process, the pressure imbalance between the chromosphere and the condensation region would push the chromospheric plasma upwards to supply mass to the corona in the form of siphon flows (Xia et al. 2011; Zhou et al. 2014). In spite of the complicated dynamics of prominence threads, it is interesting to perform the stability analysis of these structures. For this reason, we developed a very idealised model that provides the basic physics to understand the evolution of a density enhancement or thread that moves along a curved concave magnetic field, a very common situation in the solar corona. The main assumption in the model is that the behaviour of the thread is essentially that of a piston that in our case is governed by gas pressure and gravity. Our model facilitates the existence of a suspended mass in the absence of magnetic dips.

We derived two equations for the equilibrium of a thread initially located at any position along the tube (Eq. (23)) and the corresponding eigenfrequency of oscillation (Eq. (41)), which revealed useful knowledge about the stability properties of the system. In spite of the simplicity of our model the equilibrium equation is very rich from a mathematical point of view and leads to significant physical interpretations. Different scenarios are possible; a density enhancement initially located at the apex is always in equilibrium in our model, but this equilibrium is either stable or unstable. We derived analytical expressions for the bifurcation values of the different parameters when a transition between a stable and unstable regime is produced. If the thread near to the apex is unstable then the system allows two additional symmetric equilibria at lower heights that are stable. Under these circumstances a thread situated at the apex of the tube moves along the magnetic field until it reaches the lower state energy associated with the new equilibrium. When the density enhancement is not initially at the apex we have in essence a similar behaviour: a single stable solution or three solutions, two of which are stable and the other unstable. In this last case, the mathematical behaviour of the system is characterised by the appearance of a S-shaped bifurcation, well known in bifurcation theory.

Our analysis reveals the presence of bifurcation points, obtained analytically when the thread is initially located at the apex. For example, when the condition given by Eq. (35) is satisfied, the system is stable, independent of the length of the magnetic field lines where the thread is suspended. Otherwise, a thread with a fixed density contrast and length can be stable when it is located along short magnetic field lines, but unstable on long magnetic field lines, namely, higher up in the magnetic arcade. This can have some implications regarding the finite height of observed solar prominences, but also concerning the stability of unbounded models with height such as those of Hood & Anzer (1990). However, the presence of magnetic dips, missing in our model, can have important consequences regarding stability.

The length of the magnetic tubes plays a key role in the equilibrium and stability of the threads. The physical interpretation of the dependence of the instability on the length of the field lines can be connected with the pressure-volume work. The apex location for threads in short flux tubes is more stable because, with the same perturbation, the gas pressure increases more drastically in the compressed region with the decreasing volume when the initial gas volume is small (and decreases more rapidly in the expanding region).

A natural extension of the model presented in this work is the inclusion of magnetic dips through the straightforward modification of the projected gravity along the field. This most likely has relevant consequences regarding stability and affects the eventual bifurcation points that may appear in the system. This work is in progress and is addressed in part II as a continuation of the present paper part I. Another future improvement to the model is the implementation of the variation of the cross-sectional area along the tube, allowing us to incorporate a truly 3D effect in the calculation of the equilibrium and the corresponding eigenfrequencies (see Luna et al. 2016b).

An et al. (1988) and Wu et al. (1990) investigated the effects of plasma injection on the formation of the Kippenhahn-Schlüter model of prominence in optimum conditions. These authors found that for high values of the plasma-β parameter (the ratio of plasma pressure to magnetic pressure) the magnetic arcade develops a magnetic dip at the centre of the structure that supports the prominence plasma. However, comparing with our study, in the low plasma-β regime (or under others injection conditions) they found that the dip is less deep and the system develops two additional plasma enhancements located at the lateral edges of the magnetic arcade. Recent works suggest that the deformation of the magnetic field lines is determined by the parameter δ (the ratio of the gravity to the magnetic pressure) (Zhou et al. 2018; Zhang et al. 2019). An et al. (1988) suggested that the steady lateral plasma accumulates because of both the injection process and because the field lines without dips do not geometrically contain the injected plasma, but Wu et al. (1990) proposed that the prominence mass is also supported by an increase in the pressure gradient. Since in this study we consider that magnetic field lines do not change owing to the presence of the dense thread we investigate in detail the results of An et al. (1988) and Wu et al. (1990) in part II of this paper. On the contrary, studies of the formation of 1D filament threads by chromospheric heating in the presence of non-adiabatic effects, such as radiative losses and thermal conduction, show that for magnetic loops without a dip, the plasma condenses but it streams along the magnetic field and disappears after falling to the footpoints (Antiochos et al. 2000; Karpen et al. 2006). Moreover, when the thread is initially in a thermal and force-balance equilibrium state but it is disturbed by a strong velocity perturbation, the prominence mass drains down to the chromosphere. Dense blobs of falling plasma have been habitually observed (Schrijver 2001; de Groof et al. 2005), therefore it seems that threads cannot be held static along vertical magnetic flux tubes in the corona. The coronal part of the tube can only slow down the falling blobs. Müller et al. (2004) proposed that the acceleration reduces because the pressure of the cooling plasma underneath the radiating blobs slows down the descent, and Oliver et al. (2014, 2016) and Martínez-Gómez et al. (2020) argue that pressure gradient is the main force that opposes the action of gravity. Our study proposes that the pressure gradient can cause the equilibrium of threads in quasi-vertical flux tubes without dips even though it has not been corroborated by observations. Besides, our model is relatively simple to study the stability of filament threads in the magnetic field without dips. Other processes such as radiation and heat conduction must be considered, which might change the stability results.

Our model also demonstrates that when the mass of the thread is high enough the position of the equilibrium is near to the footpoints. This allows us to constrain the parameters to find solutions that do not represent suspended threads but correspond to threads that travel down the tube until they essentially settle down into a new equilibrium around the base of the corona. This last situation is applicable, for example, to coronal rain falling along magnetic field lines. In this regard, our study can be relevant to explain some results of 2D and 2.5D numerical simulations of dense blobs released in the solar atmosphere (e.g., Mackay & Galsgaard 2001; Kohutova & Verwichte 2017a,b). These simulations have, in common with our work, a dense blob (akin to our thread) that is released in a gravitationally stratified background atmosphere with a closed bottom boundary at the base of the corona. This set-up allows the dense material falling along the magnetic tube to act as a piston. For this reason, Mackay & Galsgaard (2001) and Kohutova & Verwichte (2017a,b) find, as we do, that when the blob density is large enough this cold material falls down to the surface because the increased pressure that builds underneath it cannot counteract the blob weight. On the contrary, blobs with smaller mass can find an equilibrium position before reaching the surface and so they oscillate around this position. In addition to the relevance of the blob density, ρ_t, we also note that the length of magnetic field lines and the blob length are important parameters that control the system dynamics.

Acknowledgments

A.A., J.T., R.O., and M.C. acknowledge the support from grant AYA2017-85465-P (MINECO/AEI/FEDER, UE), to the Conselleria d’Innovació, Recerca i Turisme del Govern Balear, and also to IAC³. A.A. acknowledges the Spanish ‘Ministerio de Economía, Industria y Competitividad’ for the ‘Ayuda para contratos predoctorales’ grant BES-2015-075040. We are grateful to Rafel Prohens from the Departament de Ciències Matemàtiques i Informàtica, Universitat de les Illes Balears (UIB) for his advise on bifurcation theory. All the numerical calculations in this work have been performed using Mathematica (Wolfram Research 2020).

References

All Figures

Fig. 1. Sketch of the basic model. Upper panel: situation in which the plasma thread, represented by the dark rectangle of length lt, is located at the tube apex. This configuration is in an equilibrium that can be stable or unstable. Lower panel: new equilibrium of the thread at a lower height. The curved magnetic field is uniform and directed along the tube axis. The gravity force is pointing downwards and is projected along the field line according to Eq. (1). The footpoints are supposed to be anchored at the base of the corona, represented by the shaded grey area.

In the text

Fig. 2. Left- (black line) and right- (red line) hand sides of Eq. (19) as a function of s/L, for three different values of lt. Top panel: system is linearly stable with one solution (see circle). Bottom panel: system relaxes to a new equilibrium at a lower height (see the two complementary solutions) and the trivial solution becomes an unstable solution. In all cases, L = 10 H and ρt = 100 ρ0; lt = L/110 (top panel), lt = ltb (middle panel), and lt = L/30 (bottom panel). The value ltb ≈ L/77.08 is obtained when the derivatives of the two curves at s = L/2 are the same (see Eq. (20)). The vertical dashed lines represent the edges of the domain of the variable s, namely s = lt/2 and s = L − lt/2.

In the text

Fig. 3. Solutions se of Eq. (19) as a function of ρt/ρ0 (top panel), lt/L (middle panel), and L/H (bottom panel). The expression lt = 0.1 H in the top and bottom panels and ρt/ρ0 = 100 in the middle panel. The solid lines correspond to stable solutions and dashed lines to unstable solutions. The thin vertical dashed lines correspond to the bifurcation values that have been calculated by solving Eqs. (19) and (20) simultaneously and coincide with the results in Eq. (32) (top panel), Eq. (33) (middle panel), and Eq. (34) (bottom panel).

In the text

Fig. 4. Left- (black line) and right- (red line) hand sides of Eq. (23) as a function of s/L, for three different situations with one, two, and three solutions. In these plots s0 = 0.4/L (where the black line crosses zero), L = 10 H, lt = L/50; ρt/ρ0 = 100 (top panel), ρt/ρ0 = ρtb/ρ0 (middle panel), and ρt/ρ0 = 200 (bottom panel). The value ρtb/ρ0 ≈ 137.28 is obtained from Eqs. (23) and (24). The vertical dashed lines represent the edges of the domain of the variable s, namely s = lt/2 and s = L − lt/2.

In the text

	Fig. 5. Same as Fig. 3 but for Eq. (23) and for different values of s0/L. L = 10 H and lt = 0.2 H (top panel); ρt/ρ0 = 100 and L = 10 H (middle panel); lt = 0.2 H and ρt/ρ0 = 137 (bottom panel). The bifurcation points (intersections between vertical and horizontal thin dashed lines) have been calculated by solving Eqs. (23) and (24) simultaneously.
In the text

Fig. 6. Solutions se of Eq. (23) as a function of s0/L for ρt/ρ0 = 200 (orange lines) and ρt/ρ0 = 400 (black lines). The expression L = 10 H and lt = 0.2 H. The curves correspond to $s_{\mathrm{e}}^{\mathcal{L}}$ and $s_{\mathrm{e}}^{\mathcal{R}}$ (solid lines) and $s_{\mathrm{e}}^{\mathcal{A}}$ (dashed line). The bifurcation points around s0/L = 0.3 and s0/L = 0.7 (intersections between vertical and horizontal orange dashed lines) have been calculated by solving Eqs. (23) and (24) simultaneously.

In the text

Fig. 7. ω2 as a function of ρt/ρ0 (top panel), lt/L (middle panel), and L/H (bottom panel) for $s_{\mathrm{e}}=s_{\mathrm{e}}^{\mathcal{A}}$ (below bifurcation points and lower branch curves, Eq. (31)) and $s_{\mathrm{e}}=s_{\mathrm{e}}^{\mathcal{L}}=s_{\mathrm{e}}^{\mathcal{R}}$ (upper branch curves; Eq. (40)). In these plots s0 = L/2. lt = 0.1 H in the top and bottom panels and ρt/ρ0 = 100 in the middle panel. The bifurcation values (thin dashed vertical lines) have been calculated using the analytical expressions given by Eqs. (32)–(34).

In the text

	Fig. 8. ω2 as function of different equilibrium parameters (Eq. (41)) for three different initial positions of the thread, s0. The expression L = 10 H and lt = 0.2 H (top panel), ρt/ρ0 = 100 and L = 10 H (middle panel), and lt = 0.2 H and ρt/ρ0 = 137 (bottom panel). The vertical lines correspond to the bifurcation points and are calculated by imposing that ω2 = 0.
In the text

Fig. 9. ω2 as a function of s0/L (Eq. (41)) for the three solutions obtained in Fig. 6. The solid lines correspond to $s_{\mathrm{e}}^{\mathcal{L}}$ (lines decreasing from left to right) and $s_{\mathrm{e}}^{\mathcal{R}}$ (lines increasing from left to right), and the dashed lines to $s_{\mathrm{e}}^{\mathcal{A}}$. The horizontal line represents ω = 0 and divides the plot into stable and unstable solutions. In this plot the left and right solutions are stable while the near apex solution is unstable. The bifurcation points (vertical dashed lines) have been calculated by solving Eqs. (23) and (24) simultaneously, and it is equivalent to impose ω2 = 0 in Eq. (41) once se is known.

In the text