Safety First: Stability and Dissipation of Line-tied Force-free Flux Tubes in Magnetized Coronae

Magnetized plasma columns and extended magnetic structures with both footpoints anchored to a surface layer are an important building block of astrophysical dissipation models. Current loops shining in X-rays during the growth of plasma instabilities are observed in the corona of the Sun and are expected to exist in highly magnetized neutron star magnetospheres and accretion disk coronae. For varying twist and system sizes, we investigate the stability of line-tied force-free flux tubes and the dissipation of twist energy during instabilities using linear analysis and time-dependent force-free electrodynamics simulations. Kink modes (m = 1) and efficient magnetic energy dissipation develop for plasma safety factors q ≲ 1, where q is the inverse of the number of magnetic field line windings per column length. Higher-order fluting modes (m > 1) can distort equilibrium flux tubes for q > 1 but induce significantly less dissipation. In our analysis, the characteristic pitch μ˜0 of flux-tube field lines determines the growth rate ( ∝μ˜03 ) and minimum wavelength of the kink instability ( ∝μ˜0−1 ). We use these scalings to determine a minimum flux tube length for the growth of the kink instability for any given μ˜0 . By drawing analogies to idealized magnetar magnetospheres with varying regimes of boundary shearing rates, we discuss the expected impact of the pitch-dependent growth rates for magnetospheric dissipation in magnetar conditions.

This work combines insights from different fields of plasma astrophysics.First, we exploit the vast literature on magnetic flux rope dynamics in the solar corona (e.g., Galsgaard & Nordlund 1997;Amari et al. 2003;Rugg et al. Gerrard et al. 2004;Török & Kliem 2005;Török et al. 2010;Gordovskyy & Browning 2011;Gordovskyy et al. 2014;Pinto et al. 2016;Ripperda et al. 2017a,b;Sauppe & Daughton 2018), where the injection of twist and helicity from surface motion produces a variety of dissipative events.Second, we use well-established constraints from laboratory plasma physics (e.g., Hazeltine & Meiss 2003;Bergerson et al. 2006;Longaretti 2008).The socalled safety factor denotes the inverse of the ratio of field line windings per column length and indicates the susceptibility of a flux tube to plasma instabilities (e.g., Goedbloed & Hagebeuk 1972;Goedbloed et al. 2019).Third, we repurpose numerical methods from the study of relativistic jets.Growth rates of perturbations to rotating equilibrium flux tubes of infinite length were derived numerically by Sobacchi et al. (2017), and we closely follow their implementation and analysis.Finally, we use the results from three-dimensional numerical models of a global magnetar magnetosphere to evaluate the astrophysical implications of our findings.Carrasco et al. (2019) and then with broader parameter ranges.Mahlmann et al. (2023) find eruption scenarios that range in onset time and dissipation for flux tubes twisted at one end by surface motions.However, their work does not explore a reliable instability criterion on the eruption of three-dimensional twisted flux tubes.We analyze the plasma safety factor as an instability criterion for coronal loops around compact astrophysical objects by drawing analogies to well-established theories.For example, Raadu (1972) and Hood & Priest (1979, 1981) studied the stabilizing effect of line-tying for magnetic columns with critical safety factors to explain the energy release of coronal mass ejections (CMEs).Solar flux tubes pinch regions of weak magnetic fields with shallow gradients at their edges.In contrast, magnetar current loops can have sharp edges on a strong background field, with negligible contributions by plasma pressure gradients.Their magnetic field lines are tied to a perfectly conducting stellar surface, while the ends of a solar flux tube extend through the photosphere and chromosphere.
Force-free electrodynamics (FFE), the vanishinginertia limit of ideal magnetohydrodynamics (MHD), is a good approximation for modeling the global dynamics of highly magnetized magnetospheric plasma (see, e.g., Gruzinov 1999;Blandford 2002;Komissarov 2004;Spitkovsky 2006;Parfrey et al. 2012;Carrasco & Reula 2017;Most & Philippov 2020;Yuan et al. 2020;Ripperda et al. 2021).FFE methods gain efficiency by disregarding the exact physics of nonideal dissipation, namely the screening of electric fields E ∥ along the magnetic field or in electrically dominated regions (E > B).However, especially with high-order numerical techniques, one can capture with good accuracy the evolution of magnetic pressure and tension as well as (non)linear interactions of plasma modes.In this work, we use FFE to model the growth of instabilities in perturbed flux tube equilibria with static line-tied boundaries.We probe an instability criterion for the onset of the kink mode and give limits on the dissipated magnetic energy for different evolution scenarios.
This paper is organized as follows.Section 2 semianalytically studies the linear growth of perturbations to force-free equilibrium flux tubes.Section 3 validates the expected instability growth (Section 3.1) and analyzes the dissipation of twist energy for various flux tube parameters (Section 3.2).Our discussion in Section 4 provides scalings and limits of the instability growth (Section 4.1), applies our findings to magnetized astrophysical coronae (Sections 4.2 and 4.3), and discusses limitations (Section 4.4).We state conclusions in Section 5 and give additional details on the instability evolution in Appendix A.

FORCE-FREE FLUX TUBE EQUILIBRIA
Formed by twisted magnetic field lines arched with footpoints frozen into a surface layer, coronal flux ropes are expected in several astrophysical systems, such as stars and magnetized accretion disks (see introductory references).Surface motions can drag along the line-tied magnetic field lines and disrupt the magnetized flux tube equilibria.To study the stability of highly magnetized twisted magnetic fields, we analyze simplified flux tube geometries embedded in uniform background magnetic fields (see Figure 1).We evaluate Maxwell's equations, for ideal electric fields with frozen-in magnetic flux and for a vanishing Lorentz force For stationary electromagnetic fields, Equation (6) becomes the equilibrium condition: For cylindrical coordinates (r, ϕ, z), magnetic field lines can move around the symmetry axis with a velocity v F ϕ = rΩ F , where Ω F is the field line angular velocity.In stationary and axisymmetric configurations, the only nonvanishing component of the electric field is E r = −rΩ F B z /c (Equation 5).The angular velocity and corresponding electric field are conserved along flux surfaces.In this geometry, the radial component of Equation ( 7) yields the generalized force-free Grad-Shafranov equation of a flux tube1 : We express the magnetic field B ϕ in the orthonormal basis.This work studies line-tied force-free equilibria with footpoints anchored to perfectly conducting boundaries and a vanishing field line angular velocity Ω F .In this limit, Equation ( 8) becomes We can define the inverse pitch parameter μ = B ϕ /B z and write Equation ( 9) as This equation determines the equilibrium magnetic fields of a static flux tube without field line rotation for radial profiles of the pitch μ(r).It is well studied throughout the literature in Newtonian ideal MHD (e.g., Goedbloed et al. 2019;Goldston 2020) and has exact solutions in some cases like uniform twists (Gold & Hoyle 1960) or certain oscillating magnetic fields (Lundquist 1950).In this work, we study field line columns similar to those twisted by footpoints moving on a stellar surface or an accretion disk.We choose simple pitch profiles to capture two main properties: the rapid decay of twist at flux tube boundaries, and a twist profile compatible with certain surface motions.In the following, we evaluate the stability of flux tubes with where ν ≥ 1.For a characteristic length scale r 0 , we choose f (r) = 0.5 × tanh[(r/r 0 − a)/b] with a ≈ 1.11 and b = 0.1.This choice induces a pitch profile with a maximum at r/r 0 = 1, vanishing for r ≫ r 0 .For the following instability analysis, we use ν ∈ {1, 2} to probe different pitch profiles in the flux tube.By integrating Equation (10) we generate equilibrium magnetic fields for boundary conditions B z (r ≫ r 0 ) = B bg , where B bg = B bg ẑ.We then use solutions to Equation (10) as background fields in an instability analysis of helical force-free MHD equilibria (Solov'ev 1967;Lyubarskii 1999;Sobacchi et al. 2017).

Instability analysis
We use the numerical method introduced by Sobacchi et al. ( 2017) to find growth rates for linear instabilities of flux tubes with inverse pitch profiles given by Equation (11).We numerically derive growth rates of linear perturbations ξ of the form ξ ∝ ξ(r) × e i(ωt+mϕ−kz) . (12) Here, ξ(r) is the complex-valued perturbation along the flux tube radius with frequency ω and wavenumbers (m, k) along the ϕ and z-directions, respectively.The vertical wavelength associated with such a perturbation is given by λ = 2π/k and its growth rate by the imaginary contribution Im(ω).In practice, we follow Sobacchi et al. (2017) and discretize a complex-valued balance Growth rates of the m = 1 mode for varying pitch and profiles of the twisted flux tube (Equation 11) as a function of the unstable mode wavelength λ.Circles indicate the numerically derived rates for ν ∈ {1, 2} and μ0 ∈ {0.5, 0.75, 1.0, 1.25, 1.5}.We indicate the critical length of safety factor q = 1 by vertical lines for different pitch parameters μ0.
Growth rates of the m = 2 (fluting) mode, as Figure 2.
equation for perturbations ξ to helical force-free equilibria (Lyubarskii 1999, Equation 20) on a one-dimensional mesh along the radial direction of the flux tube 2 .We impose radial boundary conditions ξ ′ (0) = 0 with an arbitrary normalization ξ(0) = 1, as well as ξ(r ≫ r 0 ) = 0 (see Section 2.2 in Sobacchi et al. 2017).We neglect line-tying boundaries of the perturbation ξ along the zdirection (see limitation in Section 4.4).For freely chosen wavenumbers (m, k), an initially estimated complex frequency ω is then driven to a solution of the balance 2 The technique of linearizing balance laws like Equation ( 8) is common in the Newtonian MHD literature (e.g., Hain & Lust 1958;Frieman & Rotenberg 1960;Goedbloed 1971).For certain background fields, growth rates of the kink mode (like Figure 2) were derived early on (e.g., Goedbloed & Hagebeuk 1972).The strategy adopted by Lyubarskii (1999) and Sobacchi et al. (2017) and employed in this work follows the analysis of nonrelativistic MHD equations for helical stationary flows by Solov'ev (1967).We refer the reader to these works for a vast background of the earlier theoretical development.
equation by minimizing the residual error of the shooting method (e.g., Vetterling & Press 1992).

The m = 1 (kink) mode
We first quantify the dynamics of the fastest-growing nonaxisymmetric instability of the flux tubes, that is, the kink mode.A commonly employed measure of the susceptibility of magnetic columns to kink instability is the so-called safety factor q.This parameter represents the inverse of the number of magnetic field line windings distributed along the tube length L: In the setup given by Equation ( 11), the safety factor is minimal at r/r 0 = 1.An instability is expected for q ≲ 1, and for each radial pitch profile we define the critical length corresponding to this threshold as L 0 = 2πr 0 /μ 0 .
In Figure 2, we display the growth rates for the kink mode in configurations given by Equation ( 11) as a function of wavelength λ m=1 = 2π/k m=1 for varying pitch profiles.We indicate the critical length L 0 for which q = 1 in Equation ( 13) by vertical lines.The actual system length then determines the instability growth.For wavelengths L < L 0 , no unstable m = 1 modes can be found.Therefore, the safety factor threshold of q = 1 is a valid criterion for the onset of the kink instability.With flux tubes long enough to allow for unstable wavelengths to fit into the system, the fastest-growing mode with L 0 < λ m=1 < L will dominate.The growth of the kink mode becomes faster for increasing pitch factor μ0 .These basic features hold for the different radial profiles of the inverse pitch with ν ∈ {1, 2}.However, the overall growth rates are significantly lower for larger ν (as a consequence of the nonlinear profile of the pitch parameter).

The m = 2 (fluting) mode
We extend the instability analysis to the m = 2 mode in Figure 3. Modes with m > 1 grow with wavelengths shorter than the critical length L 0 , below which the m = 1 (kink) mode is suppressed.Again, their maximum growth rates increase for larger values of μ0 .As in the case of the m = 1 (kink) mode, growth rates are significantly reduced for paraboloidal (ν = 2) pitch profiles when compared to linear (ν = 1) pitch profiles.Comparing the growth rates between the m = 1 mode (Figure 2) and the m = 2 mode (Figure 3), one finds ratios of Im(ω m=2 )/Im(ω m=1 ) = 0.75 − 0.93 for ν = 1 and Im(ω m=2 )/Im(ω m=1 ) = 0.58 − 0.72 for ν = 2.Although these measurements confirm the m = 1 mode as the fastest-growing instability, the growth rate of the m = 2 mode can become comparable.We evaluate the possibility of mode mixing in the following sections on simulated instability dynamics.

SIMULATIONS
We use FFE simulations to examine the instability growth in flux tubes described by Equation (11).For this, we employ a high-order FFE method with optimized hyperbolic/parabolic cleaning parameters (Mahlmann et al. 2020a,b;Mahlmann & Aloy 2021) that benefits from the Carpet driver (Goodale et al. 2003;Schnetter et al. 2004) and the Einstein Toolkit (Löffler et al. 2012;Zlochower et al. 2022) 3 .FFE simulations integrate Maxwell's equations (1-4) with currents set by the force-free condition (6) and the constraints There are notable differences between FFE and the nonrelativistic limit of MHD commonly used in laboratory and solar plasma physics.Only the drift velocity of frozen-in field lines is available in FFE; the plasma pressure and flow velocity are ordered out.Variations of the electric field can create local charge densities.Violations of the force-free constraints (14-15) rapidly dissipate nonideal electric fields.In FFE, the plasma modes reduce to Alfvénic and fast waves, both propagating with the speed of light and their characteristic polarizations (see, e.g., Komissarov 2002).In this section, we track the displacement of magnetic field lines during various instabilities and quantify the induced dissipation.
The simulations fill a rectangular domain of size x × y × z = [−4r 0 , 4r 0 ] × [−4r 0 , 4r 0 ] × [0, L] with resolution ∆x = ∆y = ∆z = r 0 /N .We choose N = 20 as the number of grid points resolving the flux tube radius.The boundaries in the xy-directions are periodic.In the z-direction we use perfectly conducting surfaces with frozen-in field lines (see Munz et al. 2000;Mahlmann et al. 2023).Simulations are initialized with background magnetic fields B z and B ϕ as solutions to Equation (10) determined by specifying the pitch profile in Equation ( 11).The employed high-order FFE method suppresses discretization noises required to seed the instability growth.Therefore, we initialize the driftvelocity perturbation where v 0 /c = 0.01, k z = 2π/L, and m = 1.In practice, this drift perturbation is set up by initializing electric fields according to Equation (5).In this configuration, the m = 1 (kink) mode dominates, and strong currents wind around the initial flux tube center.Figure 10 shows the time evolution of this setup.We first isolate the instability dynamics of force-free flux tubes for the m = 1 (kink) and m = 2 (fluting) modes.To this end, we follow the evolution of a setup with ν = 1 and μ0 = 0.75 for two different flux tube lengths a 0 ≡ L/r 0 = 8 and L/r 0 = 14.As shown in Figure 2, the kink mode of this configuration does not grow below L/r 0 ≈ 8.4, and the maximum growth rate of the m = 1 mode is captured for L/r 0 ≳ 9.2.Both configurations allow for the m = 2 (fluting) mode to grow, as shown in Figure 3.We use the projection Λ of the conserved force-free current along the magnetic field to visualize flux tubes in arbitrary field line geometries.This component of the current can be written in the form j ∥ = λB, with ∇Λ • B = 0; hence, Λ is constant along magnetic field lines.

Instability dynamics
Figure 4 shows the field line configuration and currents after the onset of instability for L/r 0 = 14 (accompanied by significant dissipation of twist energy; see Section 3.2).The setup develops clear features of the m = 1 (kink) mode, namely asymmetric variations of currents along the toroidal direction.The flux tube current cross section in panel c of Figure 4 exhibits typical structures of the kink instability (as discussed by Davelaar et al. 2020;Mahlmann et al. 2023).
The development of the kink mode is suppressed for the L/r 0 = 8 setup shown in Figure 5.The instability develops differently from the longer flux tube discussed above.Strong currents develop with an m = 2 symmetry along the toroidal direction.The characteristic fluting manifests as a thinning of the flux tube in the y-direction with bulging along the x-direction.We note that both setups evaluated in this section (L = 8 and L = 14) become unstable at similar times.As established in Section 2.1, the maximum growth rates of the m = 1 and m = 2 modes are comparable.We find Im(ω m=2 )/Im(ω m=1 ) ≈ 0.048/0.052= 0.92 (see Fig- ures 2 and 3).We study the dependence of twist dissipation on the flux tube length and simultaneous growth of the kink and higher m modes in the following section.

Dissipation of twist energy
We scan the parameter space of pitch profiles by varying ν and μ0 and measuring the dissipated magnetic energy e diss .We define the dissipated energy as the difference between twist energy before and after the development of instability, where the twist energy is We evolve perturbed equilibrium states in time for a duration of ct/r 0 = 300 − 500, making sure that the dynamical phase of the instability is fully captured and dissipation has returned to the low level of numerical diffusion.During the instability, twist energy is lost in steep gradients via numerical diffusion or by removing nonideal field components violating conditions (14-15).The total amount of dissipated energy is e diss .Figure 6 displays the twist energy dissipation as a function of the safety factor q for a set of 58 different flux tubes.For large safety factors q ≫ 1, no notable instability occurs, and dissipation is limited to numerical diffusion regardless of the pitch profile ν (empty circles).Sufficiently low safety factors q ≳ 1 allow for the growth of fluting (m = 2) or higher-order modes.However, the dissipation in this region of the parameter space remains low with e diss /e twist ≲ 0.2.For q ≲ 1 the m = 1 (kink) mode can grow, as was shown in Section 2.1.The dissipation of twist energy jumps to larger values at q ≈ 1, ranging between e diss /e twist ≈ 0.6 for ν = 1 and e diss /e twist ≈ 0.4 for ν = 2.For very low .Energy dissipated during the evolution of perturbed flux tube equilibria for various pitch profiles with ν ∈ {1, 2} and μ0 ∈ {0.5, 0.75, 1.0, 1.25, 1.5}.We display measurements of the dissipated energy e diss /etwist and a moving average (gray line).Filled circles denote setups that develop instability; empty circles represent configurations that did not show any dissipation above the numerical diffusivity.The dashed vertical lines indicate the critical safety factor (q = 1).The kink (m = 1) mode will grow for setups located to the left of this line.
safety factors q ≪ 1 the fraction of dissipated energy e diss /e twist increases further to e diss /e twist ≈ 0.8.
Configurations with q ≲ 1 develop both m = 1 (kink) and m = 2 (fluting) or higher-order modes.If the scale of the maximum growth rate of the kink instability is captured, the m = 1 mode dominates in these cases.The growth of the instabilities quickly drives the system to an event with rapid dissipation and a rearrangement into a relaxed state of lower energy.Once the fluting instability develops for q > 1, the system rearranges and dissipates energy equally fast.For q ≈ 1, when the system size allows for the development of the m = 1 mode but does not yet capture its maximum growth rate, the instability dynamics is more complex.First, m = 2 or higher-order modes develop, driving the fluting of the flux tube and mild dissipation of twist energy.At later times, the m = 1 (kink) instability significantly reduces the twist energy.Such events at the threshold of the 10 -4 0.001 0.010 0.100 Maximum growth rate of the m = 1 (kink) mode (bottom panel) and corresponding wavelength (top panel) for an extended range of pitch parameters μ0.Configurations with a very small initial twist μ0 ≪ 1 still show a growth of the kink mode.However, their maximum growth rate decays fast with Im(ωmax)r0/c ∝ μ3 0 (dashed gray line, bottom panel) for μ0 ≪ 1, and the required system length for the m = 1 (kink) mode scales with λmax/r0 ∝ 1/μ0 (dashed gray line, top panel).
critical safety factor with the subsequent development of modes of lower order can last three to five times longer than the cases of q < 1 and q > 1.

Scales and limits of the instability growth
We confirm in Section 2.1 that the m = 1 (kink) mode is the fastest-growing instability of the system that also dissipates twist energy most efficiently (Section 3.2). Figure 7 (bottom panel) extends the growth rate analysis shown in Section 2.1 for the m = 1 (kink) mode to μ0 = B ϕ /B z ≪ 1.All probed configurations, even with a small initial twist, have unstable solutions with a maximum growth rate Here, η is a parameter that depends on the radial pitch profile; it is obtained empirically as η ≈ 0.15 for ν = 1 and η ≈ 0.03 for ν = 2.We can estimate the time scale for the growth of instabilities as t i ≡ 1/Im(ω max ): By combining Equations ( 18) and ( 19) we can estimate the growth rate of the kink instability close to the critical safety factor q ≈ 1, where the flux tube aspect ratio is a 0 ≡ L 0 /r 0 ≈ 2π/μ 0 : The top panel of Figure 7 confirms the relation a 0 ∝ 1/μ 0 .The fastest-growing wavelength of the m = 1 (kink) mode requires flux tubes with aspect ratios of up to a 0 > 100 for μ0 ≪ 1.In other words, flux tubes have to be long and 'skinny' to become unstable for small values of μ0 .
The instability analysis and simulations presented in this work have magnetic equilibria (Equation 8) as a starting point.However, the time scale of instability growth has to be compared to other characteristic time scales of the system.In realistic scenarios, one relevant scale is the rate at which a flux tube is twisted by footpoint motions in the line-tied boundaries.We denote the twisting time scale as t twist and identify two relevant regimes: The slow-twisting regime with t i ≪ t twist , and the fast-twisting regime with t twist ≲ t i .In the slow twisting limit, the pitch parameter μ0 of a flux bundle of fixed length L will gradually increase until the system is disrupted by an m = 1 or higher-order instability.As the configuration slowly approaches q ≳ 1, the dissipation during instability in this regime is likely low (see Section 3.2 and Figure 6).In the fast-twisting regime, the inverse pitch μ0 can increase beyond the critical value for a fixed system size L. The safety factor can, thus, reach q ≲ 1 and significant dissipation of twist energy will likely occur (see Figure 6).Regardless of the safety factor, currents remain in the domain after instability and relaxation to a steady state.As we discuss in Appendix A, the total currents of disrupted flux tubes drop less than their total twist magnetic fields.We interpret our findings in the context of different astrophysical environments in the following sections.

Kink (m = 1) instability in the magnetar corona
Twisted magnetic fields likely play a key role in the region of active plasma processes in the relativistic magnetar magnetosphere, the so-called magnetar corona (see, e.g., Beloborodov & Thompson 2007;Beloborodov 2013;Chen & Beloborodov 2017).If their twist grows beyond a critical angle, the sheared dipole magnetosphere can erupt in flaring events with large-scale reconnection regions and energy dissipation (Parfrey et al. 2013;Mahlmann et al. 2019;Yuan et al. 2020;Sharma et al. 2023).Three-dimensional flux tubes with twisting footpoints in a disk-like patch on the magnetar surface allow for rich eruption dynamics with lateral (torus-like) and helical (kink-like) instabilities (Carrasco et al. 2019;Mahlmann et al. 2023).Figure 8 adapts a visualization of the threshold between large-scale global eruptions of the magnetosphere and confined eruptions by Mahlmann et al. (2023, red dashed line).To connect the global context of the magnetar magnetosphere to the findings of this work, we display the magnetic pitch μ0 for critical flux tubes (q = 1) of different aspect ratios.The aspect ratio a 0 = L 0 /r 0 is inferred from the length of the center field line in a dipolar flux tube as well as the radius r 0 of the disk that induces foot point motion on the stellar surface (see Figure 1 in Mahlmann et al. 2023).Flux tubes with low levels of inverse pitch require small diameters with large aspect ratios to become critical, a 0 = L 0 /r 0 = 2π/μ 0 for q = 1.Shorter dipolar flux tubes located farther away from the poles need larger diameters and pitch parameters to become critical.For the simulation parameters scanned in Mahlmann et al. (2023, gray dots in Firgure 8), instabilities occur after ct/L ≈ 25 for flux tubes closer to the poles (θ c = 30 • ) and ct/L ≈ 100 for those closer to the equator (θ c = 60 • ).In combination with the angular dependence of the radial magnetic field in the dipole magnetosphere, B r = 2µ * cos θ/r 3 , this difference in time before eruption suggests a scaling in critical pitch that is consistent4 with the contours displayed in Figure 8.
Time scales of field line displacement on the magnetar surface are not yet well constrained.They range from slow quasi-steady shearing with t twist on the order of years (ω s > 1 rad yr −1 , cf.Parfrey et al. 2012) to rapid crust deformations with millisecond creep times during flares (Thompson et al. 2017;Thompson 2022).The slow shearing is clearly separated from the instability growth time t i .However, for aspect ratios of a 0 ≳ 6.3, Equation ( 20) projects growth times above millisecond duration (also depending on the parameters η and r 0 ).Thus, rapid crust deformations can reach the fast twisting regime with t twist ≲ t i .In this limit, the safety factor can reach q < 1 due to the rapidly driven growth of pitch μ0 .As a consequence, significant part of the magnetospheric twist energy can dissipate (see Figure 6).We note that the simulations in Mahlmann et al. (2023) use ω s < 1/25 × c/R * , equating to t twist > 5.2 × 10 −3 (R * /10km) s.For large aspect ratios, or low pitch parameters μ0 , this choice of twist time scale allows for t twist ≲ t i and possibly enhances magnetospheric dissipation.We acknowledge the limitation of such direct comparison between the straight flux tubes discussed in this work and the dipolar magnetosphere in Section 4.4.

Mixed instabilities in magnetized coronae
In the case of flux bundles twisting in the slow limit of t i ≪ t twist or if the twist ceases while q > 1, m = 2 (fluting) and higher-order modes can develop.As we demonstrate in Section 2.1, the m = 2 mode grows slower than the m = 1 mode, though their growth rates are in general comparable.In the dynamic phase of the instability, this coincidence of growth rates manifests by mixing of the symmetric, short-wavelength fluting and the asymmetric, long-wavelength kinking patterns (see Appendix A).We find systems that are only susceptible to the m = 2 mode to dissipate a comparably small fraction of the twist energy (Figure 6) and maintain significant currents after relaxation (Appendix A).
However, the fluting develops with its short wavelengths and can potentially drive turbulence even when the kink mode is present.While the m = 1 (kink) mode develops predominantly around the center of the flux tube, the m = 2 (fluting) mode drives the dynamics at resonant surfaces in the outer layers.
The possibility of mode mixing in flux tubes was discussed, though not observed, in the context of the solar corona by Quinn & Simitev (2022).Instead of perturbing a equilibrium configuration with a specific profile that seeds the growth of an m = 1 (kink) mode, Quinn & Simitev (2022) rely on noise introduced by the continuous twisting of a flux tube to drive the instability.With this strategy, which is somewhat closer to the realistic flux tube evolution, they allow for the development of m > 1 modes that are usually not addressed in the literature. 5Quinn & Simitev (2022) consider the role of nonideal effects for the flux tube dynamics in MHD and find that the cumulative ohmic heating is mainly driven by the kink instability.The analysis we present in Figure 6 equally suggests that most dissipation occurs during the development of the m = 1 (kink) mode.
In contrast to Quinn & Simitev (2022) we do not find a significantly delayed onset of the kink mode due to the growth of m = 2 (fluting) patterns.This may be due to the absence of resistive/viscous effects in FFE, which effectively propagates all perturbations at the fastest velocity, the speed of light.Still, for configurations erupting at q ≈ 1 we find a notable interaction between the m = 1 and m = 2 modes (see Section 2.1).With the growth of the kink mode suppressed initially, such systems exclusively develop m = 2 (fluting) dynamics.Once the configuration is no longer in equilibrium, the flux tube changes its properties, and the m = 1 (kink) mode can grow regardless of its suppression in the initial state.Thus, for line-tied relativistic flux tubes of critical length with q ≳ 1 the m = 2 mode can drive the system to rapid dissipation by the kink instability.Some of the configurations studied in this work develop higherorder, vortex-like structures in the nonlinear phase of the instability (see, e.g., Figures 10 and 11).Such 'fingers' also appear in other studies of the relativistic kink instability (e.g., Bromberg et al. 2019;Davelaar et al. 2020); their dispersion and role for dissipation in the nonlinear phase should be evaluated further.

Limitations
The presented simulations of flux tube dynamics use the force-free limit of ideal MHD.Such FFE models are not suitable to capture the conversion of magnetic energy consistently due to the absence of information about plasma inertial properties like particle number densities or nonideal electric fields (see Mahlmann et al. 2020b;Mahlmann & Aloy 2021).The measurements of dissipation shown in Section 3.2, especially its time scales, should therefore be understood as the energy loss in the limit of rapid cooling of nonideal fields.How a magnetosphere with line-tied shear fills with plasma and how active plasma processes can dissipate the injected twist even without the onset of large-scale instabilities was previously studied in reduced dimensionality (e.g., Beloborodov & Thompson 2007).Its consistent plasma dynamics for realistic coronal geometries have to be further evaluated by particle-in-cell methods.
After a disruption by an instability, especially those occurring close to the critical value of q ≈ 1, some twist energy and currents remain in the system (see Figure 6 and Appendix A).A continuing twist of the flux tube footpoints could drive further eruptions of configurations with larger and larger twist energies (as modeled by Mahlmann et al. 2023).However, in such secondary and later events, the flux tubes are no longer in an axisymmetric equilibrium as considered in this work.Calculating the safety factor and corresponding instability criteria in nonaxisymmetric states is less straightforward.It will require careful consideration of the flux tube geometry, as in Stefanou et al. (2022) who derive generalized force-free Grad-Shafranov equilibria of magnetospheres with nontrivial twisted flux ropes.We defer studying the dynamic (in)stability of such configurations to future work.
The instability analysis presented in Section 2.1 can be improved, particularly regarding the neglect of linetying on the perturbation ξ.Contrary to typical fluxtube analyses in the solar corona, we do not impose any z-dependent boundary conditions on ξ.Line-tying with ξ vanishing smoothly at field line footpoints (e.g., Hood & Priest 1979, 1981, for the photosphere) could change the derived instability threshold and its dependence on the system size (cf.Goedbloed & Halberstadt 1994).Specifically, an eigenvalue analysis including zdependent boundaries may yield different outcomes than the purely radial one-dimensional balance equation discussed in Section 2.1.The effects of z-dependent boundaries on the stability of specific flux-tube geometries will be explored in future work.The full 3D simulations discussed in Section 3 use perfect conductor boundaries, an instantaneous line-tying.We find that the dynamics in these simulations are in good agreement with the predictions of unstable wavelengths derived by the instability analysis (Section 2.1).For small line-tying scale heights, the system length is well defined, and the onedimensional analysis in Section 2.1 provides good insights, for example, estimates of the minimum system length or the maximum instability growth rate for a given pitch factor.To produce more complete dispersion relations, future instability analyses in FFE could adopt innovative frameworks like Legolas (Claes et al. 2020).Such solvers systematically analyze the eigensystem of linearized MHD and can be extended to capture nonideal effects (De Jonghe et al. 2022;Jonghe & Kuczyński 2023).
Finally, this work focuses on cylindrical flux tubes and disregards any curvature effects experienced by bent structures commonly observed in the solar corona and expected, for instance, around magnetars and magnetized accretion flows.We only use the presented findings to qualitatively supplement models that take into account coronal geometries.MHD modes were analyzed in geometries relevant to laboratory plasmas early on, most notably in extensions to toroidal geometries (Goedbloed 1975).In parallel to the vast progress in tokamak applications, MHD models with coronal geometries built up a track record in the solar physics community (e.g., Amari et al. 2003;Gerrard et al. 2004;Török & Kliem 2005;Török et al. 2010;Gordovskyy & Browning 2011;Gordovskyy et al. 2014;Pinto et al. 2016;Ripperda et al. 2017a,b;Sauppe & Daughton 2018).There are only a few comparable works for global magnetospheric instabilities around compact objects (e.g., Carrasco et al. 2019;Mahlmann et al. 2023;Most & Quataert 2023) that have to be studied in greater detail in the future.

CONCLUSION
This paper examines the behavior of line-tied flux tubes with an axial twist of μ = B ϕ /B z in the forcefree regime.We examine perturbations of force-free flux tube equilibria both analytically (Section 2.1) and in FFE simulations (Section 3.1).Depending on the flux tube parameters, our analysis predicts and shows the development of the kink (m = 1) and/or fluting (m = 2) instabilities.To test the flux tubes' susceptibility to the m = 1 mode, we apply a stability indicator called the safety factor q = 2πr 0 μ/L, which represents the inverse of the number of azimuthal windings of the magnetic field along the flux tube length.Resulting analyses of growth timescale and energy dissipation are then ap-plied to astrophysical systems such as magnetars and the lower solar corona.
The safety factor is the key criterion to determine if line-tied FFE flux tubes without initial field line motion become kink unstable (see Figure 2).We find rapid growth of the kink instability for q < 1.For q > 1, fluting (m = 2) and higher m modes may develop when the flux tube is long enough to accommodate these modes.For q ≳ 1, an initial deformation of the flux tube by fluting modes can catalyze the development of the kink instability (Section 3.1).
We tested two flux tube pitch profiles μ(r) ∝ μ0 r ν in this work, where ν ∈ {1, 2}.Theoretically, we find the ν = 1 case generates larger growth rates for both the kink and fluting instabilities.The maximum growth rate of the fluting mode is within 60-90% of the kink mode (Figures 2 and 3), resulting in instabilities that evolve on similar timescales.However, the fluting mode dissipates only about 20% of the initial twist energy, while the kink mode dissipates between 40% and 80% (Figure 6).The maximum growth rate of the kink mode, Im(ω max )r 0 /c, scales as η μ3 0 , where η is found experimentally to be 0.03-0.15depending on the pitch profile (Figure 7).The wavelength corresponding to this fastest-growing m = 1 (kink) mode, λ max /r 0 , is proportional to μ−1 0 .The explosive release of magnetic energy in magnetospheres, including around magnetized compact objects like magnetars and black hole accretion systems, can be driven by the instability of twisted magnetic flux bundles.Similar to CMEs of the Sun, twisted flux tubes anchored to a neutron star or accretion disk can erupt and dissipate magnetic energy when their twist exceeds a critical value.In nature, the twist in a flux tube is likely established by footpoint motions of magnetic field lines in the line-tied boundary.In this paper, we develop an intuition for the onset of instabilities in line-tied force-free flux tubes and the dissipation associated with such events.We suggest that fast footpoint shearing at the line-tied boundary can tap the regime of efficient dissipation during m = 1 (kink) instabilities (q < 1).However, if the shear builds up slowly compared to the growth of the kink mode, it is likely that higher-order instabilities distort the flux tube with only little magnetospheric dissipation.This by itself can drive flux tubes to states with localized q ≲ 1 regions though the safety factor is not straightforwardly obtained in nonaxisymmetric configurations.Our simulations and growth rate analysis confirm that flux tubes of any pitch value B z /B ϕ can become kink unstable given sufficient length such that q ≲ 1 (see Section 4).The idealized analysis and simulation of isolated line-tied magnetic flux bundles presented in this paper provide an intuition for the dynamics of flux tubes and the corresponding dissipation limits that can be used in the complex modeling of radiative plasma processes around magnetars and magnetized accretion disks.As discussed in Section 2.1, the length and pitch profile of flux tubes determines the dominantly growing modes and sets the amount of dissipated twist energy e diss /e twist (Section 3.2).Short flux tubes with a safety factor q > 1 suppress the kink mode (see Figure 2).higher-order modes, like the m = 2 (fluting) mode, can still grow in such systems.The middle panel of Figure 9 shows a clear characteristic of the symmetric m = 2 (fluting) mode for a selected configuration with ν = 1, μ0 = 0.75, and L = 8.Specifically, the flux tube develops a bulging in the xz-plane and thinning in the yz-plane.Increasing the length of the same pitch configuration such that q ≲ 1 allows the growth of the m = 1 (kink) mode, as shown in Figure 10.The middle panel of Figure 10 shows the symmetric, short-wavelength bulging and thinning of the m = 2 (fluting) mode, as well as the asymmetric, long-wavelength runaway pattern of the m = 1 (kink) mode.This characteristic shape of the m = 1 (kink) mode becomes more prominent for configurations with larger initial pitch parameter μ0 , like the ν = 1, μ0 = 1.25, L = 10 setup shown in Figure 11.
With the dissipation of twist energy during the development of instabilities (see Figure 6), the considered systems relax to a new steady state.In general, this relaxed state is no longer axiymmetric or solved by the force balance in Equation ( 8).As shown in the right panels of Figures 9 to 11, currents remain in the final states of the evolution.Configurations with suppressed m = 1 (kink) modes maintain significant and ordered currents in the relaxed state (Figure 9).However, systems subject to the m = 1 (kink) instability end up with lower levels of current that extend further than the initial flux tube boundary with patches of stronger currents (Figures 10 and 11).We quantify the loss of current during the instability by measuring the change in total current λ tot = λ dV in the domain.Figure 12 shows the relative difference of initial and final currents ∆λ tot /λ tot,0 as a function of dissipated twist energy e diss /e twist for all models in Section 3.2.These measurements demonstrate that currents persist for all configurations in the final state of evolution.The change in total current is below the change of total twist magnetic fields, with ∆λ tot /λ tot,0 < e diss /e twist .
Figure 10.Currents during the evolution of a selected flux tube (ν = 1, μ0 = 0.75, L/r0 = 14).The middle panel shows structures that are characteristic of the asymmetric m = 1 (kink) mode as well as the m = 2 (fluting) mode.The m = 2 mode has shorter wavelengths than the m = 1 mode, as discussed in Section 2.1.In the relaxed state (right panel) some currents remain though currents are weaker than in the dynamic phase of shorter configurations (cf. Figure 9).We provide an animated version of this figure as supplementary material (Rugg et al. 2023b).It shows the evolution during times ct/r0 = 0 to 300; the real time duration of the animation is 5 s.Reactions of the flux tube to the instability are perceivable at ct/r0 = 80 (real time: 1 s), the most dynamic phase is at ct/r0 = 115 (real time: 2 s).
Figure 11.Currents during the evolution of a selected flux tube (ν = 1, μ0 = 1.25, L/r0 = 10).The middle panel shows structures that are characteristic of the dominating m = 1 (kink).Some currents remain in the relaxed state (right panel).We provide an animated version of this figure as supplementary material (Rugg et al. 2023c).It shows the evolution during times ct/r0 = 0 to 300; the real time duration of the animation is 5 s.Reactions of the flux tube to the instability are perceivable at ct/r0 = 60 (real time: 1 s), the most dynamic phase is at ct/r0 = 70 (real time: 1 s).

Figure 1 .
Figure 1.Schematic illustration of the configurations studied in this paper (panel b) and their astrophysical context (panel a).Coronal flux ropes, as observed on the Sun and expected around magnetars and magnetized accretion disks, are usually bent (panel a) with magnetic field lines frozen to a surface layer.We mimic such surface layers by perfectly conducting surfaces as boundary conditions in our simulations.Studying the stability of bent flux tubes on arbitrary background magnetic fields is not straightforward.We, therefore, analyze simplified flux tubes as twisted magnetic field lines embedded in a uniform background magnetic field (panel b).

Figure 4 .
Figure 4. Three-dimensional visualization of the instability of a selected flux tube (ν = 1, μ0 = 0.75, L/r0 = 14).We show views of the field-aligned current Λ along the x-axis (a), the y-axis (b), the z-axis (c), and a diagonal view (d).In this configuration, the m = 1 (kink) mode dominates, and strong currents wind around the initial flux tube center.Figure 10 shows the time evolution of this setup.
Figure6.Energy dissipated during the evolution of perturbed flux tube equilibria for various pitch profiles with ν ∈ {1, 2} and μ0 ∈ {0.5, 0.75, 1.0, 1.25, 1.5}.We display measurements of the dissipated energy e diss /etwist and a moving average (gray line).Filled circles denote setups that develop instability; empty circles represent configurations that did not show any dissipation above the numerical diffusivity.The dashed vertical lines indicate the critical safety factor (q = 1).The kink (m = 1) mode will grow for setups located to the left of this line.

Figure 8 .
Figure8.Estimate of a possible pitch parameter (μ0) distribution for critical flux tubes (q = 1) in a dipolar magnetar magnetosphere (adapted fromMahlmann et al. 2023, with gray dots denoting the parameter space explored by their global magnetospheric simulations).In a dipole field, flux tubes are parametrized by the center footpoint colatitude on the stellar surface θc, and the angular extent of the twisting region 2θT .During the instability, the dipole magnetosphere can either open up in a large-scale eruption, or energy is dissipated locally (transition roughly at the dashed red line, seeMahlmann et al. 2023).Critical flux tubes with μ0 ≪ 1 require small r0 and larger length L, and thus, footpoints closer to the poles.

Figure 9 .
Figure 9. Currents during the evolution of a selected flux tube (ν = 1, μ0 = 0.75, L/r0 = 8).The middle panel shows structures that are characteristic of the m = 2 (fluting) mode.Significant currents remain in the relaxed state (right panel).We provide an animated version of this figure as supplementary material (Rugg et al. 2023a).It shows the evolution during times ct/r0 = 0 to 300; the real time duration of the animation is 5 s.Reactions of the flux tube to the instability are perceivable at ct/r0 = 75 (real time: 1 s), the most dynamic phase is at ct/r0 = 150 (real time: 3 s).

APPENDIXA.
INSTABILITY EVOLUTION AND RELAXED STATESIt is instructive to follow the dynamics of the conserved force-free current λ throughout the instability development.