Introduction

Magnetospheric dissipation likely drives at least some of the abundant X-ray activity observed from compact object magnetospheres with active coronae, such as magnetars (e.g., Göğüş et al. 1999Göğüş et al. , 2000Göğüş et al. , 2001;;Rea et al. 2009;Rea & Esposito 2011;Kaspi & Beloborodov 2017;Esposito et al. 2020) and magnetized black hole accretion disks (e.g., Haardt et al. 1994;Di Matteo et al. 1999;Chartas et al. 2009;Uttley et al. 2014;Wilkins & Gallo 2015). The stability and dynamics of flux bundles in (highly) magnetized environments are well studied and observed, including for applications to astrophysical jets (e.g., Lyubarskii 1999;Giannios & Spruit 2006;Lapenta et al. 2006;Narayan et al. 2009;Alves et al. 2018;Bromberg et al. 2019;Davelaar et al. 2020) and the solar corona (e.g., Raadu 1972;Hood & Priest 1979, 1981;Linton et al. 1998;Lapenta et al. 2006;Kumar & Cho 2014;Florido-Llinas et al. 2020;Xu et al. 2020;Quinn & Simitev 2022). Yet, insights from the considered systems, often periodic or highly constrained by the specific astrophysical scenario, cannot answer the most fundamental questions for highly magnetized flux tubes with field line footpoints frozen (line-tied) to a stellar or disk surface boundary at both ends: If, when, and how does a flux tube become unstable, and what is the amount of dissipated magnetic energy during its instability?

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;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. 2017aRipperda et al. , 2017b;;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 so-called 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 vanishing-inertia 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 linetied 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 the Appendix.

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 references in Section 1). 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

, 5

and for a vanishing Lorentz force

6

For stationary electromagnetic fields, Equation (6) becomes the equilibrium condition:

For cylindrical coordinates (r, f, z), magnetic field lines can move around the symmetry axis with a velocity v Ff = 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 tube:foot_0

We express the magnetic field B f 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 m = f 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 ˜( ) m 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

r r a b 0.5 tanh 0 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 { } n Î 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 z B bg 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

t m kz i Here, ξ(r) is the complex-valued perturbation along the fluxtube radius with frequency ω and wavenumbers (m, k) along the f 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 ( ) w Im . In practice, we follow Sobacchi et al. (2017) and discretize a complex-valued balance equation for perturbations ξ to helical force-free equilibria (Lyubarskii 1999, Equation (20)) on a onedimensional mesh along the radial direction of the flux tube. 6We impose radial boundary conditions ( ) x¢= 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 z-direction (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 equation by minimizing the residual error of the shooting method (e.g., Vetterling & Press 1992).

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 pm = L r 2 0 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 m0 . These basic features hold for the different radial profiles of the inverse pitch with { } n Î 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 m0 . 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 (

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. 2021a(Mahlmann et al. , 2021b;;Mahlmann & Aloy 2022) 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).foot_3 FFE simulations integrate Maxwell's Equations (1)-( 4) with currents set by the force-free condition (Equation ( 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 (Equations ( 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 ´´=

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 f 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 drift-velocity 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).

Instability Dynamics

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 m = 0.75 0 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.

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

(see Figures 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 m0 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

. 17

twist bg 2

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 in Equations ( 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 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 higherorder 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 critical safety factor with the subsequent development of modes of lower order can last 3 to 5 times longer than the cases of q < 1 and q > 1.

Discussion

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

. 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 timescale for the growth of instabilities as

r c r 0.33 10 Im 1km s. 19 i 5 max 0 0 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 pm º» a L r 2 0 0 0 0 We display measurements of the dissipated energy e diss /e twist a moving average (shaded gray area). 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. . In other words, flux tubes have to be long and "skinny" to become unstable for small values of m0 .

The instability analysis and simulations presented in this work have magnetic equilibria (Equation ( 8)) as a starting point. However, the timescale of instability growth has to be compared to other characteristic timescales 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 timescale 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 m0 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 fasttwisting regime, the inverse pitch m0 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 the Appendix, 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 m0 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 twisting disk that induces footpoint 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, pm == a L r 2 0 0 0 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 Figure 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, mq = B r 2 cos r 3 , this difference in time before eruption suggests a scaling in critical pitch that is consistent 8 with the contours displayed in Figure 8. Timescales of field line displacement on the magnetar surface are not yet well constrained. They range from slow, quasisteady shearing with t twist on the order of years (ω s > 1 rad yr -1 ; see Parfrey et al. 2012) to rapid crust deformations with millisecond creep times during flares (Thompson et al. 2017;Thompson 2023). 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 m0 . As a consequence, a significant part of the magnetospheric twist energy can dissipate (see Figure 6). We note that the simulations in Mahlmann et al. (2023) use

. For large aspect ratios or low pitch parameters m0 , this choice of twist timescale 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 injection 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, shortwavelength fluting and the asymmetric, long-wavelength kinking patterns (see the Appendix). 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 (see the Appendix). 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 an 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. 9 Quinn & 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

, roughly according to the difference in contours enclosing the gray dots in Figure 8. 9 We note that the numerical work conducted for this paper initially followed a similar approach. We especially considered the late-stage evolution of flux tubes where the continuous motion of footpoints was turned off after an initial twisting episode. In such setups, we found kink-like and fluting-like dynamics, seeded purely by the boundary and developing different levels of dissipation. The long evolution times until instability onset made the systematic study of such setups prohibitive for a large parameter space.

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 higher-order, vortex-like structures in the nonlinear phase of the instability (see, e.g., Appendix 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. 2021b;Mahlmann & Aloy 2022). The measurements of dissipation shown in Section 3.2, especially its timescales, 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 were 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 the Appendix). 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. (2023) 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 line-tying on the perturbation ξ. Contrary to typical flux-tube 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 (see Goedbloed & Halberstadt 1994). Specifically, an eigenvalue analysis including z-dependent 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 one-dimensional 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; 8. Estimate of a possible pitch parameter ( m0 ) distribution for critical flux tubes (q = 1) in a dipolar magnetar magnetosphere (adapted from Mahlmann 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; see Mahlmann et al. 2023). Critical flux tubes with  m 1 0 require small r 0 and larger length L and thus, footpoints closer to the poles. Török et al. 2010;Gordovskyy & Browning 2011;Gordovskyy et al. 2014;Pinto et al. 2016;Ripperda et al. 2017aRipperda et al. , 2017b;;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 m = f B B z in the force-free 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 pm = q r L 2 0

, 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 applied 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 ˜( ) m m µ n r r 0 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, ( ) w r c Im max 0 , scales as hm 0 3 , where η is found experimentally to be 0.03-0.15 depending on the pitch profile (Figure 7). The wavelength corresponding to this fastest-growing m = 1 (kink) mode, l r max 0 , is proportional to m - 0 1 . 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 f 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. 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 ll D < e e tot tot,0 diss twist . Figure 10. Currents during the evolution of a selected flux tube (ν = 1, m = 0.75 0 (An animation of this figure is available.)

Figure 12. Change in total current Δλ tot /λ tot,0 against dissipated energy e diss /e twist for the models described in Section 3.2.