© ESO, 2011

1. Introduction

With the improvement of the observational capabilities in recent years with high spatial and temporal resolution, jets and flows are clearly observed in the solar atmosphere, such as in photospheric/chromospheric magnetic structures (Katsukawa et al. 2007; Shibata et al. 2007; de Pontieu et al. 2007; Nishizuka et al. 2008), coronal loops (Winebarger et al. 2002; Ofman & Wang 2008), prominences (Lin et al. 2003, 2005; Okamoto et al. 2007; Lin et al. 2009), or coronal holes (Cirtain et al. 2007; Scullion et al. 2009).

On the other hand, there is evidence of twisted magnetic fields in the solar atmosphere. Photospheric motions may stretch and twist anchored magnetic fields, which may lead to the consequent changes in topology in higher regions. The observed rotation of sunspots may lead to twisting of the magnetic tubes above active regions (Brown et al. 2003; Yan & Qu 2007; Kazachenko et al. 2009). Also newly emerged magnetic tubes are supposed to be twisted during the rising phase through the convection zone (Emonet & Moreno-Insertis 1996; Moreno-Insertis & Emonet 1996; Emonet & Moreno-Insertis 1998; Archontis et al. 2004; Murray & Hood 2008; Hood et al. 2009). Likewhise, some chromospheric surges show evidence of twisting fields and movement (Gu et al. 1994; Canfield et al. 1996; Jibben & Canfield 2004; Liu 2008; Liu et al. 2009). Therefore, solar magnetic tubes may have been twisted at photospheric, chromospheric, and coronal levels. Solar prominences are also supposed to be formed in a twisted magnetic field (Priest et al. 1989). However, not all works on prominence formation and equilibrium assume twisted fields, so the issue is still not clear today.

A flow-aligned magnetic field may stabilise sub-Alfvénic flows against the classical Kelvin-Helmholtz instability, while a transverse magnetic field seems to have no effect on the instability (Chandrasekhar 1961; Ferrari et al. 1981; Cohn 1983). Therefore, the magnetic field topology is crucial for the threshold of flow instability; namely, the twist of magnetic tubes may affect the instability properties of axial mass flows. Twisted magnetic tubes are also subject to the kink instability when the twist exceeds a critical value (Lundquist 1951; Hood & Priest 1979). In fact, Anzer (1968a,b) showed that all force-free fields without line-tying are unstable (regardless of whether they are smooth or have interfaces), but in the solar atmosphere the assumption of force-free fields is not always valid, so it remains important to investigate the instabilities in these models.

The magnetohydrodynamic modes of a twisted magnetic flux tube without flow can be used to study the threshold of the kink instability (Dungey & Loughhead 1954; Roberts 1956; Bennett et al. 1999; Carter & Erdélyi 2008; Soler et al. 2010). The effect of a monolithic flow along the axis has been considered by Zaqarashvili et al. (2010), who conclude that its effect is a correction on Lundquist’s stability criterion, although it can trigger the Kelvin-Helmholtz instability if the wavelength is large enough. However, this axial flow simply drags the field lines along the axis (so its main effect is to Doppler-shift the frequencies), while in the observations it is often suggested that the flows are field aligned. Our aim is therefore to study this type of twisted equilibrium configurations with plasma flowing along the field lines and to obtain their stability regimes.

2. Equilibrium configurations

We consider a magnetic flux tube with radius a and density ρ_l embedded in a uniform field environment with density ρ_e. The magnetic field inside the tube is helicoidal, while outside the magnetic field is uniform and directed along the z-axis. Both ρ_l and ρ_e are supposed to be homogeneous. The cylindrical coordinate system (r,φ,z) is used. No mass flow is present outside the tube, so the surrounding coronal medium is considered to be uniformly magnetised (B_eẑ is a constant), uniform (ρ_e constant), and lacking mass flow at the equilibrium.

To study this equilibrium we use the set of ideal MHD equations, since no dissipative effects are taken into account. These equations are (1)plus the solenoidal condition ∇·B = 0. We have labelled the gas pressure as p_g to avoid confusion with other parameters that appear later in this paper.

Next we describe the equilibrium quantities (labelled 0). The equilibrium magnetic field inside the tube has the form (2)with B₀ = B_l, representing a twisted flux tube. We keep B₀ in this section to point out that these relations can be applied to other configurations. Assuming therefore that the equilibrium gas pressure at the boundary is the same inside and outside, the continuity of the total pressure across the tube boundary gives (3)This choice of equilibrium magnetic field assumes that there is a current sheet at the loop boundary.

Moreover, the plasma is assumed to flow along the field lines inside the magnetic tube, so we prescribe the equilibrium velocity field as (4)with v₀ = v_l for our model with a twisted field inside the tube, and again we keep v₀ in this section. Notice that v₀  ×  B₀ = 0, ∇·B₀ = 0, and ∇·v₀ = 0, satisfying the induction equation, the solenoidal condition and the mass continuity equation, respectively. However, the equation of motion, (5)indicates that the magnetic field is not force-free inside the tube, since there is an electric current (6)so the Lorentz force must be balanced by the advection and pressure terms. Assuming a stationary equilibrium and using some vector identities, we cast Eq. (5) as (7)This equation can also be written to introduce the total pressure, p_T0 = B₀·B₀/(2μ) + p_g0, (8)Since B₀ and v₀ are parallel vectors, this means that the pressure terms must be balanced by a tension force that is lessened by a factor  (with the Alfvén speed squared defined as ) because of the field-aligned flow.

Finally, using our equilibrium magnetic field and velocity (Eqs. (2) and (4)), this yields a differential equation for the radial functions α(r) and γ(r), (9)This gives a relation between both r-dependent components of the magnetic field and the equilibrium gas pressure: two are freely chosen, and Eq. (9) then gives the other dependence so that the equilibrium relations are satisfied.

In this paper we are interested in studying a tube with uniform pressure, so the equilibrium pressure does not depend on the radial coordinate, so p_g0 = p_gl = p_ge  =  const. Then, Eq. (9) gives a relation between both components of the equilibrium magnetic field, so we just can freely choose one of them.

2.1. Uniform axial field

A particular case of the previous section is to assume that γ(r) = 1, so the axial components of both the equilibrium velocity and magnetic field are uniform. Under these conditions, Eq. (9) gives (10)The non-trivial solution is (with K_a an arbitrary constant) (11)Following Eq. (11), there is a family of solutions whose radial dependence is prescribed by the relation between the magnitudes of the equilibrium velocity and magnetic field. Sub-Alfvénic flows (v₀ < c_A0) are only in balance if the radial dependence diverges when r → 0, but super-Alfvénic flows (v₀ > c_A0) lead to finite values at the tube axis. The particular case leads to a linear profile (uniform twist of the flux tube and j₀ = 2pB₀/μ   e_z), similar to the one considered by Zaqarashvili et al. (2010), but with the plasma flowing along the field lines instead of a purely axially monolithic flow. This is also an extension of the classical study in Chandrasekhar (1961), which is performed assuming the plasma is incompressible.

3. Wave equations

3.1. Wave equation inside the magnetic tube for uniform twist

The equations governing the dynamics of the plasma for small perturbations (labelled 1) from the equilibrium state are obtained from the linearisation of Eqs. (1), namely (12)(13)(14)(15)c_s being the sound speed. We also consider a situation in which the plasma beta is low, and then the gradient of the gas pressure perturbation in Eq. (12) can be neglected in front of the magnetic and advection terms, so Eq. (15) is not needed and we are left with seven unknown functions, namely the components of the perturbed magnetic field and the perturbed velocity and the perturbed density.

The system of equations in Eqs. (12)−(15) are a generalisation for field-aligned flows of the well-known equations for twisted tubes, which are recovered if we set v₀ = 0 (Hain & Lüst 1958; Appert et al. 1974). The equation for the perturbed density (Eq. (14)) needs to be considered with the rest, since the perturbed density now appears in the equation of motion (Eq. (12)).

Next, we assume a simple dependence for the equilibrium, α(r) = r/p, which satisfies Eq. (11) if K_a = 1/p, with p the pitch of the magnetic field, so 2πp is the longitudinal distance between consecutive turns in the twisted flux tube. Then, (16)so the flow must be super-Alfvénic to satisfy the motion equation. This equilibrium magnetic field, B₀ = B_l(0,r/p,1), is the same as the one in Zaqarashvili et al. (2010), but here the flow is field-aligned instead of being axis-oriented. The equilibrium condition across the tube surface is then .

As the unperturbed parameters only depend on the r coordinate, the perturbations can be Fourier-analysed with exp [i(mφ + kz − ωt)] . It is convenient to non-dimensionalise the perturbed quantities as (17)Then, the wave equations can be cast as (24)Next we non-dimensionalise the variables as (25)with the frequency including the factor v_lk.

We can reduce the system to two first-order differential equations after some algebra, This gives the following set of differential equations for the two remaining unknown radial functions b_z and v_r: (31)where  {A_ij}  are algebraic real coefficients depending on the parameters P, m and the radial coordinate d (see the Appendix). Both variables have a phase shift of π/4 for stable modes because of the i factor in Eqs. (31).

These two first-order differential equations can be transformed into one second-order linear differential equation for v_r with non-constant coefficients, (32)(33)The expressions of these coefficients in terms of the equilibrium parameters are given in the Appendix. Both coefficients are quotients of polynomials of high order in s; hence, the resulting differential equation is very difficult to solve analytically, and numerical methods may prove useful to obtain its solutions. In fact, Eq. (33) is similar to the Hain-Lüst equation (Hain & Lüst 1958) for twisted tubes without flow, but instead Eq. (33) is valid in a tube filled with uniform twist, low-beta plasma and a field aligned flow.

It is interesting to check the singularities of this equation, since they provide information about the solution. In addition to a singularity as s → ∞, Eq. (32) is singular if (34)which are the zeroes of the denominator of the coefficient C₂ (Eq. (A.6)). The solutions of Eq. (34) give three frequencies near which infinite modes are clustered (also called bands), (35)(36)and the position of additional singular points to infinity and 0, (37)or (38)These additional singular points only appear if Ω > Ω_accum2 or Ω < Ω_accum3 and for high values of s outside the tube, so they are not relevant for our analysis.

The position of the first band (Eq. (35)) can be related to the location of mode’s rational surfaces, which for a typical force-free field are found when the axial wavenumber satisfies (39)with no dependence on position in our equilibrium choice. Then, ω_accum1 = (1 + m/(kp))   kv_l represents this accumulation point about which the band is centred. The second and third bands (Eq. (36)) are similar to this one, but related to the Alfvén frequency being close to the Doppler-shifted frequency, since in our choice of equilibrium with uniform twist , so ω_accum2,3 =  [1 + m/(kp)]    k   (v_l ± c_Al). These bands are similar to those found in Zaqarashvili et al. (2010).

3.2. Solution outside the flux tube

So far, we have derived a differential equation for the perturbed velocity inside the twisted flux tube (Eq. (32)), but to obtain solutions we need to consider the boundary conditions at the tube boundary r = a.

First of all, the solution outside the flux tube is relatively simple, since the outside plasma is uniform and permeated by a constant magnetic field B_e. In these conditions a Bessel equation can be obtained for the total pressure perturbation (Edwin & Roberts 1983; Cally 1986; Díaz et al. 2002) after Fourier-analysing all the variables except r, (40)(), with the radial perturbed velocity related with the total pressure perturbation as (41)The solution for r > a can be written as (42)with and A_e and arbitrary constant. In our non-dimensional units (43)(44)(45)There is another important point to consider: the presence of the external magnetic field B_e has the appearance of a cut-off frequency as a consequence, since λ_e must be real to have an evanescent solution for r → ∞. This condition gives the value of the cut-off frequency (46)Finally, there is an interesting limit that is relevant to some solar atmospheric structures: high-density ratio χ. In coronal loops this ratio is above unity, but can be much higher in other structures, for example in prominences of typically χ ≈ 100. We thus consider χ ≫ 1, so the external Alfvén speed c_Ae in Eq. (44) increases. Then λ_e/k ≈ 1, the cut-off frequency (Eq. (46)) becomes high, and Eq. (43) can be simplified to (47)which for A ≪ 1 gives (48)

3.3. Boundary conditions (uniform twist)

Finally, we need boundary conditions to match these outer solutions to the inner ones. The first one is the continuity of the normal component of the Lagrangian displacement, and this condition is not modified by the presence of twist or flows. The second condition is the continuity of the total pressure perturbation for untwisted tubes, but for twisted tubes there are extra terms involving the equilibrium magnetic field azimuthal component (Dungey & Loughhead 1954; Roberts 1956; Bennett et al. 1999). However, in our model there are other extra terms coming from the presence of non-uniform flows. The conditions must be imposed on the Lagrangian total pressure, which is (Goedbloed 1983; Goedbloed & Poedts 2004; Goedbloed et al. 2010) (49)We use the convention of denoting  [a]  = a₂ − a₁ as the jump in the quantity a across the boundary, the subscripts 1 and 2 indicating the values on either side. Our boundary is r = a, so  and this equation becomes (50)Using the equilibrium condition (Eq. (8)), we obtain (51)which with our choice of constant z-component in the equilibrium magnetic field and uniform twist, becomes (52)with a plus sign in the term of the displacement.

We need to relate the displacement to the perturbed velocity, which is not trivial in a flowing plasma. Following Goossens et al. (1992); Goedbloed & Poedts (2004); Goedbloed et al. (2010), after neglecting non-linear terms the relation is (53)which using Eq. (4) gives (54)The quantity ω_d = v₀αm/r + v₀kγ = v₀(1 − m/p) is the Doppler shift in frequency. Hence, in our variables the radial components of the perturbed velocity and the displacement are related as (55)In terms of our functions, the total pressure perturbation inside the tube can be written as (56)The expressions of the coefficients E₁ and E₂ are given in the Appendix.

As a result, we have the solution outside the flux tube and the matching conditions, so we can numerically solve the differential equation. However, we also need to impose another condition: regularity at r → 0, which implies that p_T(r = 0) = 0. Then we can integrate from r = r₀ to r = a, where r₀ is a small number, chosen to avoid numerical instabilities at r = 0. Next we apply the boundary conditions in Eq. (52), to obtain the value of the constant A_e in Eq. (42), and then we use the second boundary condition to obtain the frequencies of the stationary modes of this model using a shooting and matching numerical scheme.

3.4. Thin tube limit

In many cases in the solar atmosphere we are dealing with thin tubes, that is, tubes whose length is much more longer than their width. Some approximations may exploit this condition to obtain analytical approximations. In our problem, the relevant parameter is A, whose order of magnitude is generally around 10^-2.

Taking this into account, we can use the Taylor expansion of the coefficients in Eq. (32) for s ≪ 1 and m ≠ 0, (57)(58)Again, the expressions of the coefficients D₁ and D₂ are given in the Appendix. Notice that D₂ still has the three singularities that lead to the accumulation frequencies and bands (Eqs. (35), (36)), but not the singular points in Eq. (37), since they only appear for large s. If we take the first-order approximation, Eq. (32) is written as (59)whose solution is (60)with K₁ and K₂ arbitrary constants. Since we expect the solutions to be regular as s → 0, we demand that K₁ = 0. Then for kink modes (| m| = 1) this gives v_r ≈ const., while for fluting modes (| m| ≥ 2) this gives v_r ≈ 0, confirming the boundary conditions in Sect. 3.3. The sausage mode (m = 0) is a special case, so is considered in the next section.

We can now use the following order in the expansions in Eqs. (57), (58), so the differential equation is cast as (61)which is a form of the Whittaker differential equation, whose solutions can be written in terms of confluent hypergeometric functions F(a,b;x). The general solution of Eq. (61) for m ≠ 0 is (62)Again, the condition of having regular solutions as s → 0 is satisfied if we impose K₃ = 0 (m > 0) or K₄ = 0 (m < 0). The solution for m > 0 is therefore (63)and the equivalent one for m < 0.

We can expand this expression further for small s to gain insight into the solution (64)Thus, we recover the first-order expansion results and obtain a correction.

Fig. 1 Spatial plots of χr (arbitrary units) inside the magnetic tube for m = 1. The solid line is the numerical solution of Eq. (32), the dotted line represents the thin tube approximation in terms of confluent hypergeometric functions (Eq. (63)), and the dashed line its first-order expansion (Eq. (64)). The upper panel is calculated with P = 0.5 and Ω = 1.5, while the lower panel is done with P = 2.0 and Ω = 4.0.

We test these approximations by plotting the numerical solution of the system, the solution in terms of the confluent hypergeometric functions (Eq. (63)) and its first-order series expansion (Eq. (64)) for a typical set of values for the parameters. This is seen in Fig. 1. For s < 0.1 both approximations are quite accurate, so the approximations can be used in the range appropriate to solar atmospheric conditions, while for higher values of s the analytical approximation is no longer valid. Also, for high values of P, the range where the analytical solution is accurate increases, but at the same time, this range decreases for high values of Ω.

Next, we need to write the expressions for the boundary conditions (Eq. (52)) in terms of v_r and its derivative (65)with the expressions of the coefficients given in the Appendix. Finally, using Eqs. (43) and (63) we can obtain a dispersion relation (66)with (67)Focusing on the kink m = 1 mode this expression gives (68)We can have an approximation to this dispersion relation expanding in series for a small A. Using the properties of the confluent hypergeometric functions and the approximation in Eq. (48) for the external solution, the first-order approximation of Eq. (66) for m ≠ 0 gives a sixth-order algebraic equation for Ω (69)whose coefficients are given in the Appendix. There is no simple form for expressing the roots of this equation, but it may provide a fast way to compute the solutions, even the ones with complex frequencies that may appear as solutions of this algebraic equation.

3.4.1. Sausage modes

We have so far investigated the modes with m ≠ 0. However, the modes that satisfy m = 0 (also called sausage modes) must be dealt with separately. We follow the same procedure by expanding the coefficients in Eq. (32) for s ≪ 1 and m = 0, (70)(71)Our differential equation is written in the first-order approximation as (72)whose solution is (73)with K₁₀ and K₂₀ as arbitrary constants. Since we expect the solutions to be regular as s → 0, we require that K₂₀ = 0. Now we have v_r ≈ K₁₀s, which is different from the m ≠ 0 modes. In fact, to obtain these solutions numerically the condition that needs to be imposed is v′(s) = const. This is somewhat expected, since even for a straight magnetic flux tube without twist or flow, the sausage modes behave differently as r → 0.

We can now use the following order in the expansions in Eqs. (70), (71), so the differential equation is cast as (74)which is a form of the Whittaker differential equation, whose solutions can be written in terms of confluent hypergeometric functions F(a,b;x). The regular solution as s → 0 is (75)We can expand this expression further for small s to gain insight into the solution (76)again, recovering the first order expansion results and obtaining a correction.

Finally, we can write the dispersion relation for the sausage modes as Eq. (66), but with (77)

4. Stability of the model

Once we have discussed the differential equations and boundary conditions necessary for solving them, we can look for their solutions. First of all, we need to consider which parameters are relevant to our study. After non-dimensionalising the equations we are left with a set of parameters  {m,P,A,χ,Ω} , since the flow speed and the equilibrium magnetic field are related in our equilibrium by . We first focus on the kink mode (m = 1), since it is the easiest mode to detect in observations and the one that leads to unstable modes in similar equilibrium configurations.

Fig. 2 Dispersion diagram (Eq. (66)) plotting the non-dimensional frequency ω/(vlk) against the non-dimensional tube radius ka for the values of the parameters shown in the legend. Stable solutions are plotted in solid lines, while unstable solutions are plotted in dotted lines. The numerical solution of the system in Eq. (31) has been overplotted as diamonds. The dashed lines represent the external cut-off frequencies (Eq. (46)).

Fig. 3 Zoom of the dispersion diagram of Fig. 2 (Eq. (31)) near the bands for the values of the parameters shown in the legend. Only the first ten modes farther from each accumulation frequency have been plotted in each panel. The upper panel is centred on the first band (Eq. (35)), the central panel is around the second band (positive sign of Eq. (36)), and the lower one is around the third band (negative sign of Eq. (36)). The values of the accumulation frequencies have been overplotted as a dashed line.

We first study the dependence of the frequency on the tube radius ka (Fig. 2) for a given set of the rest of parameters. The thin-tube approximation shows the right behaviour even for ka > 0.1, although the negative frequencies are not obtained correctly. There is a rich array of modes, and we can distinguish the three bands whose frequencies are clustered around the accumulation frequencies in Eqs. (35), (36), which for these parameters correspond to Ω_accum1 = 1.63, Ω_accum2 = 2.77, and Ω_accum3 = 0.48. Only a few modes in each band have been represented. A zoom of the dispersion diagram around the bands is shown in Fig. 3, where we can see that the first band is symmetric with respect to Ω_accum1, while the second and third bands only have modes on one side of their accumulation frequency and are much more clustered. The spatial shape of these modes, shown in Fig. 4, contains an increasing number of wiggles near the tube boundary as the mode frequency approaches the accumulation frequency. This can be understood better by checking previous results: in a twisted tube without flow (Bennett et al. 1999), the solution can be expressed in terms of Bessel functions, and the band appears when the argument becomes large, allowing many extrema in the interior of the tube boundary. In our model we can see a similar behaviour in the thin-tube limit: for the hypergeometric functions in Eq. (63) the D₂ coefficient becomes large as Ω → Ω_accum (Eq. (A.10)), so many extrema are added to the eigenfunction inside the tube as we approach the accumulation frequency (Fig. 4).

We can also see an instability being triggered in Fig. 2: two modes become complex at the bifurcation point at ka = 0.42, ω/(v_lk) = 1.79. These modes are the first ones in each band, with a single maximum in r = 0 and no minimum.

Fig. 4 Spatial plots of vr (arbitrary units) for the values in the legend. Only the five modes whose frequency is farthest from the accumulation frequency have been represented (Fig. 3, upper panel), each mode with a different linestyle. The vertical dashed line is the position of the tube boundary.

The dependence on the parameters can be sumarised as follows: modifying the density ratio χ affects the frequencies slightly, but raising it also increases the cut-off frequency (Eq. (46)), allowing more modes to become trapped. Even one or more bands might become leaky for low values of χ. Raising the tube radius ka expands the bands where infinite modes are clustered and eventually can lead to destabilisations, while changing the pitch kp modifies the position of the bands.

We next concentrate on the unstable modes (Fig. 5). As ka is decreased these modes merge, and beyond the bifurcation point they become unstable: now their frequencies have a complex part. For this particular choice of parameters the model would therefore be unstable for ka < A_crit = 0.152. In accordance with Zaqarashvili et al. (2010), the instability threshold is no longer located in Ω = 0 (as in Bennett et al. 1999) because of the presence of a flow, but in this case the exact value of the threshold is calculated numerically. Solving the comparatively simple Eq. (69) for the thin tube gives quite accurate results.

Fig. 5 Dispersion diagram plotting ω/(vlk) against kp for the values of the parameters shown in the legend. The upper panel represents the real part of the frequencies and the lower panel the imaginary part. Stable solutions are plotted in solid lines, while unstable solutions are plotted in dotted lines. The dashed line represents the external cut-off frequency (Eq. (46)), while the dot-dashed lines represent the position of the accumulation frequencies of each band (Eqs. (35), (36)). Only the modes that go through the bifurcations have been plotted. The results of the algebraic equation in Eq. (69) have been overplotted as dots.

Fig. 6 Stability diagram (χ = 15) for different values of the azimuthal number m. The dashed line represents the threshold without flow (Dungey & Loughhead 1954), while the dotted line represents the threshold with a monolithic axial flow with (Zaqarashvili et al. 2010). Equilibrium configurations whose parameters are above the stability lines are unstable.

Finally, we represent the threshold of the instability (Fig. 6). We can see that the stability region is very limited, and for modes with kp < 1 the tube is unstable unless ka ≪ 1. We can compare it with the results of a twisted tube without flow (Dungey & Loughhead 1954; Roberts 1956; Bennett et al. 1999) and with a twisted tube with monolithic axial flow (Zaqarashvili et al. 2010) to check that the stability threshold on the corrected Lundquist’s instability criterion, which in our variables is (78)and this type of flow lowers it substantially. The m = 1 mode is the first mode that becomes unstable as ka is increased.

5. Conclusions

We have studied the equilibrium properties of twisted magnetic flux tubes with flows along the field lines and found the relations between the equilibrium quantities that must be satisfied. One of the main differences with equilibriums studied in previous works is that the pressure must be balanced by a lessened tension force due to the advection term in the equation of motion, but otherwise the equilibriums are similar to those of twisted flux tubes without flow (see for example Priest 1982).

We then focused on one particular case, a tube with uniform twist across the radial coordinate, B₀ = B_l(0,r/p,1), since it is a particularly interesting case that has been considered in many studies. We found that the flow velocity must be super-Alfvénic to satisfy the equilibrium relations, . Although these high values can be matched with a few observations, it is interesting to point out that is the only equilibrium in a zero-beta plasma that has a uniform twist and a field-aligned plasma velocity.

Next we studied the stability of this equilibrium by obtaining the normal modes, both numerically by solving the system of ordinary differential equations obtained from the linearised MHD system and obtaining analytical approximations in the thin tube limit, which led to confluent hypergeometric functions. One needs to be careful deriving the appropriate boundary conditions in this situation, since the continuity of the perturbed velocity and the total pressure are no longer valid. The analytical approximation is accurate enough if the tube radius is small. Compared with the case of monolithic axial flow (Zaqarashvili et al. 2010), there are three bands of infinite modes clustered around the Doppler-shifted frequency with the flow speed and the flow plus and minus the Alfvén speed. This appears even with our configuration of the flow, which does not lead to a simple frequency Doppler shift in the differential equations.

Finally, we investigated the stability of the equilibrium and found that the stability region is significantly lowered. This contradicts the results in Chandrasekhar (1961) for the incompressible case. If we compute the ratio between the magnetic energy and the total energy (Eq. (297), Chapter XII) we obtain (79)which is below 1/2, and thus stable in the incompressible range following their results, but unstable in the compressible range if ka is not very small. Another important difference with previous results is that the marginal stability does not occur when ω = 0, as in Dungey & Loughhead (1954); Roberts (1956); Bennett et al. (1999), or at the Doppler-shifted frequency, as in Zaqarashvili et al. (2010). This kink instability of the model may be caused by the large magnitude of the flow needed to assure that the uniform twisted tube is in equilibrium. It is therefore necessary to relax some assumptions in the equilibrium to study equilibriums with sub-Alfvénic flows, either by having a radial dependence for the magnetic field or taking the plasma pressure into account. The investigation of these cases is left for later work.

This instability in super-Alfvénic twisted flows may play a relevant role in chromospheric surges. This plasma eruptions may reach propagating speeds of 50−200 km s^-1, over the typical Alfvén speed in the lower layers of the atmosphere. Moreover, some surges fall back or are ejected as CMEs, but sometimes they are disrupted (Jibben & Canfield 2004; Liu 2008). For example, with the data in Jibben & Canfield (2004) for NOAAA AR 7117 surge in 1992 March 28, the axial velocity is v_l = 150 km s^-1, while a rotational velocity of Ω_rot ≈ 10^-3 radian s^-1 leads to a pitch of p = v_l/Ω_rot ≈ 1.5  ×  10⁵ km. We can estimate from the panels that the surge width is about 20′′, so we obtain a radius of a ≈ 14   000 km, which gives a/p = A/P ≈ 0.1. This line has been overplotted in the stability diagram in Fig. 7. Thus, this surge would be unstable under perturbations with P < 1.7, that is, for k ≈ 10^-5 km^-1 and smaller wavenumbers. Unfortunately, it is difficult to estimate the axial wavelength of the perturbation observationally, but it is possible that this kind of instability may cause the mentioned disruption, along with this surge and other observed ones. In Liu (2008) it is suggested that the so classified diffuse surges are spread because of the overlying magnetic configuration in which the surge is propagating, but this disruption could also be caused by some instability. Again, we need to refine our analytical work to allow further comparisons with observations, but the different stability regimes could explain the behaviour of super-Alfvénic surges and the disruption of some observed surges.

Fig. 7 Stability diagram for m = 1, with the unstable region filled with dots. The dashed line represents the equilibrium models that satisfy the observational parameters in Jibben & Canfield (2004).

Finally, it is worth remarking that the unstable regions are larger in this model than in previous ones, mainly because the high value of the magnetic Mach number needed in the equilibrium to have a field-aligned flow with a uniform longitudinal speed. However, it is also known in twisted force-free fields that photospheric line-tying is an essential ingredient to stabilise the loops (Raadu 1972; Hood & Priest 1979). Other important ingredients are the current layer at the interface, which in this work has been addressed in the boundary conditions, here, other effects might be important, such as dissipation by resistivity and resonant absorption. However, these effects are expected to have little influence on the stability of the modes. Also, other equilibrium configurations need to be addressed: for example, Eq. (11) for other profiles of B_0ϕ and different values of the magnetic Mach number, or Eq. (9) for tubes without uniform longitudinal component B_0z. The stability analysis of these modes is left for later work on this topic.

Acknowledgments

A.J.D. thanks the Spanish MICINN for support under a Juan de la Cierva Postdoc Grant. A.J.D., R.O., and J.L.B. also acknowledge the financial support from the Spanish MICINN, FEDER funds, under Grant No. AYA2006-07637. R.S. acknowledges support from a Marie Curie Intra-European Fellowship within the European Commission 7th Framework Program (PIEF-GA-2010-274716). The authors also thank J. Andries for useful suggestions and the unknown referee for his/her comments and remarks, which have definitely helped us to improve the paper, and for spotting a mathematical error.

References

Online material

Appendix A: Coefficients in the equations

The coefficients that appear in Eq. (31) can be computed with the help of a symbolic computation program. These coefficients turn out to be rather cumbersome quotients of polynomic expressions in our model with uniform twist, (A.1)(A.2)(A.3)(A.4)Next, we can compute the coefficients in Eq. (33), (A.5)(A.6)Then, the coefficients that appear in the boundary conditions (Eq. (56)) are given by (A.7)(A.8)The coefficients that appear in the thin tube limit for the differential equations (Eqs. (57), (58)) can be cast as (A.9)(A.10)

The coefficients for the boundary conditions suitable for analytical solution (Eq. 65), which also appear in the dispersion relation (Eq. 66), are (A.11)(A.12)Finally, the coefficients that appear on the sixth-order algebraic equation (Eq. 69) for m > 0 are

(A.13)

All Figures

Fig. 1 Spatial plots of χr (arbitrary units) inside the magnetic tube for m = 1. The solid line is the numerical solution of Eq. (32), the dotted line represents the thin tube approximation in terms of confluent hypergeometric functions (Eq. (63)), and the dashed line its first-order expansion (Eq. (64)). The upper panel is calculated with P = 0.5 and Ω = 1.5, while the lower panel is done with P = 2.0 and Ω = 4.0.

In the text

Fig. 2 Dispersion diagram (Eq. (66)) plotting the non-dimensional frequency ω/(vlk) against the non-dimensional tube radius ka for the values of the parameters shown in the legend. Stable solutions are plotted in solid lines, while unstable solutions are plotted in dotted lines. The numerical solution of the system in Eq. (31) has been overplotted as diamonds. The dashed lines represent the external cut-off frequencies (Eq. (46)).

In the text

Fig. 3 Zoom of the dispersion diagram of Fig. 2 (Eq. (31)) near the bands for the values of the parameters shown in the legend. Only the first ten modes farther from each accumulation frequency have been plotted in each panel. The upper panel is centred on the first band (Eq. (35)), the central panel is around the second band (positive sign of Eq. (36)), and the lower one is around the third band (negative sign of Eq. (36)). The values of the accumulation frequencies have been overplotted as a dashed line.

In the text

	Fig. 4 Spatial plots of vr (arbitrary units) for the values in the legend. Only the five modes whose frequency is farthest from the accumulation frequency have been represented (Fig. 3, upper panel), each mode with a different linestyle. The vertical dashed line is the position of the tube boundary.
In the text

Fig. 5 Dispersion diagram plotting ω/(vlk) against kp for the values of the parameters shown in the legend. The upper panel represents the real part of the frequencies and the lower panel the imaginary part. Stable solutions are plotted in solid lines, while unstable solutions are plotted in dotted lines. The dashed line represents the external cut-off frequency (Eq. (46)), while the dot-dashed lines represent the position of the accumulation frequencies of each band (Eqs. (35), (36)). Only the modes that go through the bifurcations have been plotted. The results of the algebraic equation in Eq. (69) have been overplotted as dots.

In the text

	Fig. 6 Stability diagram (χ = 15) for different values of the azimuthal number m. The dashed line represents the threshold without flow (Dungey & Loughhead 1954), while the dotted line represents the threshold with a monolithic axial flow with (Zaqarashvili et al. 2010). Equilibrium configurations whose parameters are above the stability lines are unstable.
In the text

	Fig. 7 Stability diagram for m = 1, with the unstable region filled with dots. The dashed line represents the equilibrium models that satisfy the observational parameters in Jibben & Canfield (2004).
In the text