The Evolution of Writhe in Kink-Unstable Flux Ropes and Erupting Filaments

The helical kink instability of a twisted magnetic flux tube has been suggested as a trigger mechanism for solar filament eruptions and coronal mass ejections (CMEs). In order to investigate if estimations of the pre-eruptive twist can be obtained from observations of writhe in such events, we quantitatively analyze the conversion of twist into writhe in the course of the instability, using numerical simulations. We consider the line tied, cylindrically symmetric Gold-Hoyle flux rope model and measure the writhe using the formulae by Berger and Prior which express the quantity as a single integral in space. We find that the amount of twist converted into writhe does not simply scale with the initial flux rope twist, but depends mainly on the growth rates of the instability eigenmodes of higher longitudinal order than the basic mode. The saturation levels of the writhe, as well as the shapes of the kinked flux ropes, are very similar for considerable ranges of initial flux rope twists, which essentially precludes estimations of pre-eruptive twist from measurements of writhe. However, our simulations suggest an upper twist limit of $\sim 6\pi$ for the majority of filaments prior to their eruption.


Introduction
The m = 1 kink mode or helical kink instability (hereafter KI) is a current-driven, ideal magnetohydrodynamic (MHD) instability. It occurs in a magnetic flux rope if the winding of the field lines about the rope axis (the twist) exceeds a critical value (e.g., Shafranov 1957;Kruskal & Tuck 1958;Freidberg 1982;Priest 1982). The instability lowers the magnetic energy of the flux rope by reducing the bending of field lines, which leads to a characteristic helical deformation (writhe) of the rope axis. Such writhing is often observed in erupting filaments or prominences in the solar corona (Figure 1), which has led to the suggestion that the KI can trigger filament eruptions and CMEs (e.g., Sakurai 1976;Sturrock et al. 2001;Török & Kliem 2005;Fan 2005).
The KI has been studied extensively for laboratory plasmas (see, e.g., Bateman 1978;Goedbloed et al. 2010, and references therein). In applications relevant to the low-β solar corona, typically force-free, cylindrically symmetric flux rope configurations of finite length are considered. The anchoring of coronal loops and prominences in the solar surface is modeled by imposing line tied boundary conditions at the flux rope ends. Properties of the KI such as the instability threshold and growth rate, as well as the formation of current sheets, have been investigated for various radial twist profiles in both straight and arched flux rope geometries (e.g., Hood & Priest 1981;Mikić et al. 1990;Baty & Heyvaerts 1996;Gerrard et al. 2001;Török et al. 2004). MHD simulations of kink-unstable flux ropes have been employed to model coronal loop heating and bright-point emission (Galsgaard & Nordlund 1997;Haynes & Arber 2007), soft X-ray sigmoids , energy release in compact flares (Gerrard & Hood 2003), microwave sources in eruptive flares , and rise profiles, rotation, and writhing of erupting filaments and CMEs (Török & Kliem 2005;Williams et al. 2005;Fan 2005; Kliem et al. 2012). In spite of this large body of work, the amount and evolution of the writhing in kink-unstable flux ropes was quantified only very rarely (Linton et al. 1998;Török et al. 2010). Systematic investigations of the dependence of the writhe on parameters such as the initial flux rope twist or geometry have not yet been undertaken.
The quantity writhe measures the net self-coiling of a space curve. It is related to the total torsion along the curve: the sum of writhe and total torsion remains constant under deformations, unless the curve develops an inflexion point, where curvature vanishes (Moffatt & Ricca 1992). Twist and writhe of a thin flux rope are related to its magnetic helicity via H = F 2 (T + W ), where F is the axial magnetic flux, T is the number of field line turns, and W is the writhe of the rope axis (Cȃlugȃreanu 1959;Berger & Prior 2006). The writhe for flux ropes with footpoints on a boundary (such as the photosphere) can be defined by the same formula, using relative helicity for H (Berger & Field 1984). W depends only on the shape of the axis of the rope; while T measures the net twist of the field lines in the rope about the axis. Since magnetic helicity is conserved in ideal MHD, the KI converts twist into an equal amount of writhe. Here we quantify this process for the first time systematically for a range of initial flux rope twists, using MHD simulations.
For our study we consider the straight, uniformly twisted, force-free flux rope equilibrium by Gold & Hoyle (1960, hereafter GH), line tied at both ends. In the absence of knowledge about typical twist profiles in coronal flux ropes and due to its force freeness, the equilibrium serves as a convenient reference model. Mechanisms other than the KI that may cause writhing (see Kliem et al. 2012 for a detailed discussion) are excluded. Furthermore, the KI of the GH model does not lead to the formation of a helical current sheet, which triggers reconnection in the nonlinear development of other flux rope equilibria (e.g., Baty & Heyvaerts 1996;Gerrard & Hood 2003). Therefore, the evolution of the axis deformation can be followed well into the saturation phase of the writhe, which makes this equilibrium particularly suited for our purpose. We measure the axis writhe using the formulae by Berger & Prior (2006), which express the quantity as a single integral in space, facilitating its calculation.
Our motivation for this study is derived from the interest in obtaining estimates of the twist in pre-eruptive solar configurations from the amount of writhing observed during an eruption. At present, the twist cannot be obtained directly, since the magnetic field cannot be measured in the coronal volume and since extrapolations from photospheric vector magnetograms are not yet sufficiently reliable in practice, especially for volumes containing a filament (McClymont et al. 1997;Schrijver et al. 2008). Twist estimations based on the observations of pre-eruptive coronal configurations are hampered with substantial uncertainties (see Section 4). The writhe of erupting filaments, on the other hand, can be obtained with a reasonable accuracy if the filament displays a coherent shape ( Figure 1) and if observations from more than one viewing angle are available (for example from the STEREO mission; Kaiser et al. 2008) or if the eruption is directed toward the observer . Although twist estimates from the writhe can only be obtained in retrospect, they may facilitate systematic studies of this possibly critical parameter for CME initiation and may be useful for comparison with other means of estimation.

Flux Rope Model and Numerical Setup
The GH model used in this study is given by and represents a uniformly twisted, force-free flux rope of infinite radial extent. The constant b is related to the axial length of one field line turn, Λ, often referred to as the pitch, by b = 2π/Λ and, at the same time, represents the inverse scale length of the radial field profile. Using the customary form of the expression for the twist angle, a GH rope of length L has a twist of Φ = bL. This is related to the number of field line turns T = L/Λ by Φ = 2πT . In our calculations we fix the scale b = π, so that the radial field profiles are identical in all simulations, and we vary the initial twist by varying the flux rope length.
The numerical set-up for the simulations is the same as in our previous studies of the kink and torus instabilities (e.g., Török & Kliem 2005; for a detailed description see Török & Kliem (2003). The compressible ideal MHD equations are integrated using the simplifying assumptions of vanishing pressure and gravity on a discretized Cartesian box with uniform spacing ∆ = 0.04 and L x = 8, L z = 16. In order to take advantage of the z-axis line symmetry inherent in the configurations considered, we orient the flux rope parallel to the y axis, so that the integration need be carried out only in the "half box" {y > 0}. The flux rope length equals the full box length 2L y . Except for {y = 0}, where mirroring according to the z-axis line symmetry is applied, the MHD variables are held fixed at their initial values at all boundaries (and for consistency the velocity is kept at zero also one grid layer inside the boundaries). This models line-tying at the ends of the rope, |y| = L y .
We vary the initial twist, Φ 0 = 2πL y , in the range (3-10.6)π by varying L y in our series of runs. Since the rope expands strongly and in different ways in the different runs ( Figure 2), we minimize the influence of the top and bottom boundaries by positioning the rope axis at appropriate initial heights z = h 0 . See Table 1 for the values of L y , Φ 0 , and h 0 used in this investigation.
The initial density distribution is specified to be ρ 0 = B 3/2 0 , such that the Alfvén velocity decreases slowly with distance from the flux rope axis (which corresponds to the conditions in the solar corona). The MHD variables are normalized by quantities derived from a characteristic length of the initial equilibria, chosen to be Λ/2, and the initial magnetic field strength and Alfvén velocity at the flux rope axis. All runs start with the fluid at rest. A small initial velocity perturbation localized at the flux rope center is imposed in all runs (analogous to Török et al. 2004).

Results
Since the chosen initial twists in the series all exceed the KI threshold for the line tied GH equilibrium of Φ cr ≈ 2.5π (Hood & Priest 1981), all configurations are unstable. The helical nature of the growing perturbation is clearly visible in the linear phase (the early phase of the instability during which the exponentially growing amplitude of the axis displacement remains small; top two rows in Figure 2). In the nonlinear phase (when the amplitude of the axis displacement becomes large) the flux rope starts to expand strongly by the action of the hoop force, which comes into play as soon as the flux rope develops some overall curvature between its line tied ends (bottom two rows in Figure 2 and Section 4).

Flux rope axis evolution
We measure the growth of writhe at the axis of the flux rope using Equations (4)-(6) in Török et al. (2010); see also Berger & Prior (2006). Note that flux surfaces away from the axis undergo a smaller deformation, with less twist converted into writhe, but Figure 2 indicates that a substantial cross section of the GH rope attains similar writhe. For our strongly twisted cases, the measurements are reliable only until the perturbed flux surfaces start to approach the boundaries of the box.
Figure 3(a) shows the development of writhe by the KI. The writhe grows exponentially in the linear phase and then reaches saturation in the nonlinear phase of the instability. It can be seen that the saturation level of the writhe does not scale linearly with the initial twist, which is different from what one might intuitively expect.
The flux ropes with the smallest twist, Φ 0 = 3π and 4.5π, exhibit a very similar behavior, except for a significantly faster initial evolution of the run with Φ 0 = 4.5π. The writhe saturates at W ≈ 0.95 in both runs, corresponding to a converted twist of Φ ≈ 1.9π in the vicinity of the flux rope axis. The morphological evolution and the resulting axis shapes are very similar too (Figure 2). In both cases, the axis deforms into a one-turn helix. Figure 3(b) shows that the evolution of writhe coincides well with the release of magnetic energy by the KI and the displacement of the flux rope axis.
A somewhat different evolution takes place for Φ 0 = 7.5π: the writhe first reaches a maximum after the initial exponential growth phase, then decreases by ≈ 20 percent, but subsequently starts to increase again slowly, reaching W ≈ 1.05 at the end of the simulation. Figure 3(b) shows that the flux rope continues to rise (at a slower rate) after the writhe has reached its first maximum, accompanied by ongoing magnetic energy release. The energy saturates when the writhe reaches its temporary minimum. Some further release occurs later on in the evolution; this appears to be related to reconnection that occurs at outer flux surfaces when those approach the boundary of the simulation box. The morphological evolution is somewhat different from the smaller-twist cases: the axis shape obtained in the nonlinear phase of the instability is still dominated by a one-turn helix, but it becomes internally helically deformed (see Figure 2). The decrease of the writhe appears to be related to the reversal of the orientation of the flux rope legs in the vicinity of the ±L y boundaries (see the second and fourth panel for Φ 0 = 7.5 π in Figure 2), which occurs at the same time as the writhe decrease. This is a consequence of the line-tying that would be absent in infinitely extended flux ropes, and it implies a temporary increase of twist in the vicinity of the flux rope axis. The subsequent increase of the writhe is most likely associated with the development of the internal axis deformation, but may be related to some degree also to the reconnection mentioned above.
The run with Φ 0 = 6.0 π is an intermediate case: both the temporary decrease of the writhe and the internal helical axis deformation are present but rather weak, and the writhe at the end of the simulation (W ≈ 1.0) lies in between the corresponding values for the smaller-twist runs and the run with Φ 0 = 7.5π.
The cases with the largest twists (Φ 0 = 9π and 10.6π) show a very different behavior. After the fast initial rise, the writhe continues to grow at a much smaller rate until it reaches a maximum of W ≈ 1.75 − 1.8, after which it slowly decreases (the decrease is not visible for Φ 0 = 9π in the figure, since the evolution during the nonlinear phase is significantly slower than for Φ 0 = 10.6π). The flux rope axis now develops a helix with about two turns (see Figure 2 for Φ 0 = 9 π). We attribute the slow increase and decrease of the writhe to the complex morphological evolution of the flux rope for large twists: the development of two expanding helices within the finite domain forces approaching flux rope sections to give way to one another -an effect that is much less pronounced in cases where only one helix develops. The decrease may be also a consequence of the strongly expanding helices approaching the boundaries of the simulation box.

What determines the amount of twist converted into writhe?
From Figure 3(a) it is obvious that the conversion of twist into writhe depends on the initial twist in a non-trivial manner. Figures 2 and 3 indicate that it is related to the number of helical turns the rope axis develops in the course of the instability. This number depends on the wavelength of the most unstable mode and on the range of unstable modes permitted by the finite length of the rope. In the linear phase of the instability, the helical eigenmode with the highest growth rate dominates the way the rope starts to deform. The finite length of the line tied rope modifies this picture in the nonlinear phase.
The growth rate as a function of axial wavelength for the GH equilibrium is shown in Figure 4 (from Linton et al. 1998). This plot is for the case of infinite axial but finite radial extent of the rope, R = 3π/2b. Linton et al. (1998) find the peak growth rate and its location to be unchanged for larger R (our x-z box sizes correspond to R ∼ 8π/b). The growth of the helical kink mode peaks at the wavelength λ = 1.85Λ, meaning that the KI grows fastest at a writhe wavelength of about twice the twist pitch. For our choice b = π we have 2Λ = 4. Therefore, for a double helix to dominate in the linear phase, the required box length is 2L y > 8, equivalent to a required twist Φ = 2πL y > 8π.
Although this corresponds nicely to the jump of the final writhe in Figure 2 between twists Φ 0 = 7.5π and 9π, it does not actually explain the occurrence of a jump. From the stability analysis we know that the most rapidly growing mode at λ ≈ 2Λ is permitted to occur as soon as the box length satisfies 2L y > 2Λ = 4. Therefore, for 2L y > 2Λ = 4 non-integer values for the number of turns of ≈ L y /Λ can occur at the peak growth rate in the linear phase. This agrees with the simulation results shown in the two upper rows of Figure 2. The dominant mode in this phase exhibits a little more than one turn for Φ 0 = 4.5π, nearly two turns for Φ 0 = 7.5π and a little more than two turns for Φ 0 = 9π.
The jump in the writhe can only be understood from the nonlinear evolution of the instability. This shows a clear tendency to develop an integer number of turns. As long as two axial wavelengths don't fit into the box, the mode with a single turn dominates in this phase. Contributions of the linearly most strongly growing mode, which then has a shorter wavelength, are clearly present (most obviously for Φ 0 = 7.5π), but no longer dominate. We attribute this result to the action of the hoop force for kinking flux ropes of finite length. The line-tying leads to an axial dependence of the displacement which is absent for infinitely extended ropes. For Φ 0 ≤ 7.5π the displacement is largest in the mid-plane {y = 0} of our symmetric simulations and tapers off toward the line tied ends at y = ±L y (see the two upper rows in Figure 2). An overall net bending of the rope in the direction of the displacement in the mid-plane (which points along the z axis in all our simulations) results. This introduces a Lorentz self-force in the rope, known as hoop force, pointing in the direction of the bending (Bateman 1978). The hoop force amplifies the perturbation of the rope axis in the mid-plane above its purely KI-driven displacement, thus creating a one-turn helix for Φ 0 ≤ 7.5π. Such amplification occurs at a pair of symmetrically located displacements for Φ 0 ≥ 9π, creating a helix with two turns. We discuss the implications of this result for filament eruptions in the following section.

Discussion
We studied the conversion of twist into writhe in a simulation series of the KI in the GH model, considering initial flux rope twists in the range 3.0 π ≤ Φ 0 ≤ 10.6 π. We found in all cases a saturation of writhe in the nonlinear phase of the instability, after an initial exponential increase during the linear phase. However, the final writhe does not scale simply with the initial flux rope twist. Rather, the amount of twist converted into writhe seems to be determined predominantly by the number of helical turns the flux rope axis develops in the nonlinear phase.
For 3.0 π Φ 0 7.5 π, the rope axis develops a one-turn helix. For twists close to the upper end of this range, internal helical deformations of the one-turn helix develop, due to helical eigenmodes with wavelengths λ < 2L y . However, the axis shape remains to be dominated by one turn in the nonlinear phase of the instability. The resulting writhe is close to unity for all cases, corresponding to a converted twist of ∼ 2 π in the vicinity of the flux rope axis. If the twist is increased beyond this range, the rope axis develops more than one helical turn, and considerably more twist is converted into writhe (∼ 3.5 π in our simulations with Φ 0 = 9.0 π and Φ 0 = 10.6 π). We attributed the relatively similar writhe values obtained in each respective range, as well as the pronounced increase of the writhe between them, to the action of the hoop force on line-tied, kink-unstable flux ropes of finite length.
The basically discontinuous dependence of the final writhe upon the initial twist displayed in Figure 3 essentially precludes a reasonable estimation of the initial twist from observations of the writhe in solar filament eruptions and CMEs. The saturation levels of the writhe are very similar for initial twists up to Φ 0 ≈ 8π, requiring an accuracy of writhe determination for such an estimate that cannot be reached in solar observations. Moreover, the final writhe, as any other property of the KI, depends on the radial twist profile of the initial equilibrium. Therefore, a precise knowledge of this profile, combined with a parametric simulation study like the one in Figures 2 and 3 for a range of different profiles, would be required to permit a reliable estimate of twist. Further effects of importance for the final writhe enter when arched flux rope equilibria are considered (see Török et al. 2010), rendering a twist estimation from writhe observations even more difficult.
In order to compare our results with the KI in force-free equilibria with non-uniform radial twist profile, we performed simulation series similar to the one presented here for the straight flux rope model termed "Equilibrium 2" in Gerrard et al. (2001) and the arched flux rope model by Titov & Démoulin (1999). Unfortunately, in all runs the flux rope axis was destroyed by reconnection at current sheets before the writhe would clearly saturate (see Amari & Luciani 1999, Haynes et al. 2008, and Valori et al. 2010 for examples of such reconnection), so that these simulations cannot be used for the purpose of this study.
However, our writhe measurements for the KI in the GH equilibrium provide at least a rough upper limit for the initial twist of erupting filaments. It is observed that kinking filaments typically do not display more than one helical turn and hardly any significant internal helical deformation of their axis. Combined with our simulations, this suggests that the initial twist typically does not exceed values Φ 0 ∼ 6π. This is supported by simulations of the KI in the Titov-Démoulin model, which show strong internal helical deformations for twists above this value (see, e.g., Figure 1 in Kliem et al. 2010 and Figure 12 in Török et al. 2010).
Occasionally, however, the Sun seems to succeed in building up higher twists. Several examples can be found in Vršnak et al. (1991), whose estimates of the end-to-end twist fall in the range (3-15) π for a sample of prominences close to the time of eruption. A particularly clear indication of very high twist (of ∼ 10π) was obtained by Romano et al. (2003) for a filament eruption on 19 July 2000 (Figure 1c). These estimations are based on measurements of the pitch angle of selected helical prominence threads, which are then converted into twist assuming a uniform radial twist profile both along and across the axis of the underlying flux rope. The latter assumption may be a severe oversimplification, since force-free flux ropes embedded in potential field must generally have a nonuniform radial twist profile in order to match the field at the surface of the rope (see, e.g., Figure 2 in Török et al. 2004). Still, the simulations presented here support the existence of such a high twist at least for the case shown in Figure 1c, based on the strong bending in the lower part of the filament legs (see also Figure 12b in Török et al. 2010, where a strongly nonuniform radial twist profile was used). A further observed case of very high twist may have been an apparently three-fold helix described in Gary & Moore (2004).
While these estimations remain uncertain to a considerable degree, we can ask how such large twists, if present, may be produced in the solar corona. It is widely believed that twist is accumulated prior to an eruption by flux emergence (e.g., Leka et al. 1996), photospheric vortex flows (e.g., Romano et al. 2005), or the slow transformation of a sheared magnetic arcade into a flux rope (e.g., Moore et al. 2001;Aulanier et al. 2010). It has been argued that flux ropes that form by one or more of these mechanisms will become kink-unstable long before the large twists mentioned above can be reached. While this is likely true for the majority of cases, several scenarios for the build-up of large twists appear to plausible. First, the KI threshold can vary in a wide range as a function of the thickness of the rope (e.g., Hood & Priest 1979;Baty 2001;Török et al. 2004), so sufficiently thin flux ropes may be able to harbor large twist in a stable state. Second, sufficiently flat highly twisted flux ropes may be stabilized by strong ambient shear fields , or by gravity if they contain sufficient filament material. Third, significant twist may be added by reconnection to the rising flux in the course of an eruption (e.g., Qiu et al. 2007). Finally, flux ropes may reconnect and merge prior to an eruption, thereby adding up their respective twists (e.g., Pevtsov et al. 1996;Canfield & Reardon 1998;Schmieder et al. 2004;van Ballegooijen 2004).
We thank the anonymous referees for very helpful suggestions and Z. Mikić for providing a routine that was helpful for the writhe calculations. The contribution of T.T. was supported by NASA's HTP, LWS, and SR&T programs and by the NSF. M.G.L. received support from NASA/LWS and the ONR 6.1 programs. B.K. was supported by the DFG. The research leading to these results has received funding from the European Commission's Seventh Framework Programme under the grant agreement No. 284461 (eHEROES project). L.vDG.'s work work was supported by the STFC Consolidated Grant ST/H00260X/1 and the Hungarian Research grant OTKA K-081421.  : Writhe (diamonds), vertical displacement of the flux rope axis from its initial position on the z axis (dash-dotted lines; scaled to fit into the plot), and total magnetic energy (dashed lines; normalized to initial value) for Φ 0 = 4.5 π (red) and Φ 0 = 7.5 π (green). Fig. 4.-Helical kink instability growth rates ω vs. axial wavelength Λ/λ for the Gold-Hoyle equilibrium of infinite axial but finite radial extent of R = 3π/2b (from Linton et al., 1998; c AAS. Reproduced with permission). b = B θ /(rB z ) is the twist per unit length and v A is the initial Alfvén velocity at the flux rope axis.