In the solar atmosphere flux ropes are subject to currentdriven instabilities that are crucial in driving plasma eruptions ejections and heating. A typical ideal magnetohydrodynamics instability developing in flux ropes is the helical kink which twists the flux rope axis. The growth of this instability can trigger magnetic reconnection which can explain the formation of chromospheric jets and spicules but its development has never been investigated in a partially ionized plasma PIP. Here we study the kink instability in PIP to understand how it develops in the solar chromosphere where it is affected by chargeneutral interactions. Partial ionization speeds up the onset of the nonlinear phase of the instability as the plasma  of the isolated plasma is smaller than the total plasma  of the bulk. The distribution of the released magnetic energy changes in fully ionized plasma and PIP with a larger increase in internal energy associated with the PIP cases. The temperature in PIP increases faster also due to heating terms from the twofluid dynamics. PIP effects trigger kink instability on shorter time scales which is reflected in more explosive chromospheric flux rope dynamics. These results are crucial to understanding the dynamics of smallscale chromospheric structuresminifilament eruptionsthat thus far have been largely neglected but could significantly contribute to chromospheric heating and jet formation.
['Solar chromosphere', 'Solar magnetic reconnection', '1479', '1504', 'Astrophysics - Solar and Stellar Astrophysics', 'Physics - Plasma Physics']
Kink Instability of Flux Ropes in Partially Ionized Plasmas
{'Kink instability of flux ropes in partially-ionised plasmas': '<!-- image -->  \nG/i.pc/u.pc/l.pc/i.pc/a.pc  \nM/u.pc/r.pc/t.pc/a.pc/s.pc  \n,  \n1, 2  \nA/n.pc/d.pc/r.pc/e.pc/w.pc  \nH/i.pc/l.pc/l.pc/i.pc/e.pc/r.pc  \n,  \n2  \n/a.pc/n.pc/d.pc  \nB/e.pc/n.pc  \nS/n.pc/o.pc/w.pc  \n2  \n1 Los Alamos National Laboratory, Los Alamos NM, 87545, USA  \n2 Department of Mathematics and Statistics, University of Exeter, Exeter EX4 4QF, UK', 'ABSTRACT': 'In the solar atmosphere, flux ropes are subject to current driven instabilities that are crucial in driving plasma eruptions, ejections and heating. A typical ideal magnetohydrodynamics (MHD) instability developing in flux ropes is the helical kink, which twists the flux rope axis. The growth of this instability can trigger magnetic reconnection, which can explain the formation of chromospheric jets and spicules, but its development has never been investigated in a partially-ionised plasma (PIP). Here we study the kink instability in PIP to understand how it develops in the solar chromosphere, where it is affected by charge-neutral interactions. Partial ionisation speeds up the onset of the non-linear phase of the instability, as the plasma 𝛽 of the isolated plasma is smaller than the total plasma 𝛽 of the bulk. The distribution of the released magnetic energy changes in fully and partially-ionised plasmas, with a larger increase of internal energy associated to the PIP cases. The temperature in PIP increases faster also due to heating terms from the two-fluid dynamics. PIP effects trigger the kink instability on shorter time scales, which is reflected in a more explosive chromospheric flux rope dynamics. These results are crucial to understand the dynamics of small-scale chromospheric structures - mini-filament eruptions - that this far have been largely neglected but could significantly contribute to chromospheric heating and jet formation.  \nKeywords: Solar chromosphere (1479) - Magnetohydrodynamics (1964) - Solar magnetic reconnection (1504)', '1. INTRODUCTION': 'In the solar chromosphere, the formation of highly dynamic features such as spicules (Roberts 1945) and jets (Mohler 1951) is strongly bound to changes in connectivity of the magnetic field. The energy fuelling such explosive events can be released through the onset of magnetohydrodynamics (MHD) instabilities, that promote the conversion of stored magnetic energy through processes such as magnetic reconnection (Yamada et al. 2010; De Pontieu et al. 2012). The energy release associated with these processes promotes the formation of several types of dynamical structures, and can provide an important contribution to the chromospheric heating (Roberts 1945; Mohler 1951; Osterbrock 1961).  \nFlux ropes, which are ubiquitous in the solar atmosphere, are twisted by convective motions, reaching an equilibrium that can be unstable to the kink modes (Hood & Priest 1979; Hood 1992; Priest & Forbes 2000). Kink modes define a class of ideal, current-driven MHD instabilities that can become particularly important in triggering explosive events. When subjected to the kink instability, a plasma column bends and develops displacements that are transverse to the column axis and driven by an unbalance in magnetic pressure (Hood 1992). The instability depends both on the length of the flux tube and the twist of magnetic field lines wrapped around it (Velli et al. 1990): if the field lines wrap multiple times along the length of the flux tube, then the system becomes unstable to kink modes (Biskamp 2000).  \nWhen a plasma column is traversed by an axial magnetic field, the 𝑚 = 1 kink mode develops a helical deformation of the plasma column. Such deformation is magnified by the compression of the magnetic field lines at the internal side of the bends, with a consequent increase in magnetic pressure, and the expansion of field lines at the external side of the bends, that results in a decrease of magnetic pressure (Goedbloed & Poedts 2004). The magnetic pressure imbalance leads to an increase in the distortion  \ngiuliamurtas31994@gmail.com  \nof the plasma column, which consequently results in a larger magnetic pressure gradient and a faster growth of the instability. The onset of the helical kink instability reduces the magnetic energy within the flux rope by converting the twist of the magnetic field lines in an helical deformation of the rope axis itself, namely the writhe (Torok et al. 2014). Frequent observations of writhing in erupting filaments and prominences in the solar corona led to the suggestion, supported by many works, that the helical kink instability could be responsible for triggering filament eruptions and CMEs (Sakurai 1976; Sturrock et al. 2001; Torok & Kliem 2005; Fan 2005).  \nIn systems where magnetic diffusivity is significant the kink instability develops in flux ropes, with an initial growth phase followed by a nonlinear resistive phase, characterised by current sheets formation and the onset of magnetic reconnection (Parker 1988; Hood et al. 2009; Bareford & Hood 2015; Snow et al. 2018). The compression and expansion of magnetic field lines triggers the relaxation process (Taylor 1974, 1986), where reconnection straightens the magnetic field lines and lead to the release of the energy stored in the highly stressed and twisted field configuration. The magnetic field tends to relax toward a state of minimum magnetic energy defined by a force-free field, called the Taylor state (Taylor 1974), where the magnetic helicity is preserved. This phase results in a less twisted configuration of the flux rope. During relaxation, magnetic energy is converted into kinetic energy and heat (Heyvaerts & Priest 1984; Browning et al. 2008). The role of relaxation has been examined through several studies involving numerical simulations of cylindrical coronal loops (Browning et al. 2008; Botha et al. 2011; Gordovskyy et al. 2016; Pinto et al. 2016).  \nFlux ropes have been well studied in the fully-ionised coronal environment, where the kink instability has been identified in several observations (Srivastava et al. 2010; Liu et al. 2016), but it is unclear how these structures evolve in the lower atmosphere. Below the corona, twisted magnetic fields reconnect in regions where the plasma composition is heterogeneous and the plasma is only partially-ionised: this could be the case of micro-filament and mini-filament eruptions (Sterling et al. 2015, 2016; Sterling & Moore 2016; Samanta et al. 2019; Sterling et al. 2020), taking place in small-scale loop systems observed in extreme ultraviolet and X-rays (Harden et al. 2021; Madjarska et al. 2022; Sterling et al. 2022). The kink instability has been investigated as an explanation for jet formation through several numerical studies between the upper chromosphere (Nishizuka et al. 2008) and the lower corona (Yokoyama & Shibata 1995; Moreno-Insertis & Galsgaard 2013; Fang et al. 2014). However, numerical models that study unstable flux ropes in partially-ionised chromospheric plasma are still absent. Therefore, it is essential to investigate how the helical kink instability grows in a partially-ionised plasma, and the differences with the development of the same instability in a fully-ionised plasma.  \nThe goal of this study is to examine the kink instability developing in single loops as a potential explanation for mini-filament eruptions in the solar chromosphere, and compare a fully-ionised plasma case (MHD) with two partially-ionised plasma cases (PIP). In Section 2, the models of our simulations are presented. In Section 3 results are discussed, with a focus on the energy release during the onset of magnetic reconnection within the nonlinear phase of the kink instability. Conclusions are presented in Section 4.', '2. METHODS': 'The simulations presented in this work are performed through the (PIP) code (Hillier et al. 2016; Snow & Hillier 2021). The equations are solved throughout the domain through a fourth-order central difference scheme, and the time updates are computed by a four-step Runge-Kutta scheme for time integration.  \nThe partially-ionised plasma is modelled through a two-fluid ion-neutral hydrogen plasma. Ions and neutrals are described by two separate sets of equations and coupled though elastic collisions, charge-exchange, collisional ionisation and recombination. The neutral fluid is governed by:  \n𝜕𝜌 𝑛 𝜕𝑡 + ∇ · ( 𝜌 𝑛 v 𝑛 ) = 𝐷, (1)  \n𝜕 𝜕𝑡 ( 𝜌 𝑛 v 𝑛 ) + ∇ · ( 𝜌 𝑛 v 𝑛 v 𝑛 + 𝑝 𝑛 I ) = R , (2)  \n𝜕𝑒 𝑛 𝜕𝑡 + ∇ · [ v 𝑛 ( 𝑒 𝑛 + 𝑝 𝑛 )] = 𝐸, (3)  \n𝑒 𝑛 = 𝑝 𝑛 𝛾 -1 + 1 2 𝜌 𝑛 𝑣 2 𝑛 , (4)  \n𝑇 𝑛 = 𝛾 𝑝 𝑛 𝜌 𝑛 , (5)  \n𝜕𝜌 𝑝 𝜕𝑡 + ∇ · ( 𝜌 𝑝 v 𝑝 ) = -𝐷, (6)  \nwhile the ionised fluid is governed by:  \n𝜕 𝜕𝑡 ( 𝜌 𝑝 v 𝑝 ) + ∇ · 𝜌 𝑝 v 𝑝 v 𝑝 + 𝑝 𝑝 I -BB + B 2 2 I ! = -R , (7)  \n𝜕 𝜕𝑡 𝑒 𝑝 + 𝐵 2 2 ! + ∇ · [ v 𝑝 ( 𝑒 𝑝 + 𝑝 𝑝 ) - ( v 𝑝 × B ) × B + 𝜂 (∇ × B ) × B ] = -𝐸 -Φ 𝐼 + 𝐴 heat , (8)  \n𝜕 B 𝜕𝑡 - ∇ × ( v 𝑝 × B -𝜂 ∇ × B ) = 0 , (9)  \n𝑒 𝑝 = 𝑝 𝑝 𝛾 -1 + 1 2 𝜌 𝑝 𝑣 2 𝑝 , (10)  \n∇ · B = 0 , (11)  \n𝑇 𝑝 = 𝛾 𝑝 𝑝 2 𝜌 𝑝 . (12)  \nThe subscripts 𝑝 and 𝑛 identify physical quantities of the ion-electron plasma and of the neutral fluid respectively. The variables v , 𝑝 , 𝜌 , 𝑇 and 𝑒 are the fluids velocity, gas pressure, density, temperature and internal energy, 𝛾 = 5 / 3 is the adiabatic index and B is the magnetic field.  \nThe terms 𝐷 , R and 𝐸 are respectively the source terms for mass, momentum and energy transfer between the two species:  \n𝐷 = Γ rec 𝜌 𝑝 -Γ ion 𝜌 𝑛 , (13)  \nR = -𝛼 𝑐 𝜌 𝑛 𝜌 𝑝 ( v 𝑛 -v 𝑝 ) + Γ rec 𝜌 𝑝 v 𝑝 -Γ ion 𝜌 𝑛 v 𝑛 , (14)  \n𝐸 = -𝛼 𝑐 𝜌 𝑛 𝜌 𝑝 " 1 2 ( v 2 𝑛 -v 2 𝑝 ) + 𝑇 𝑛 -𝑇 𝑝 𝛾 ( 𝛾 -1 ) # + 1 2 ( Γ rec 𝜌 𝑝 v 2 𝑝 -Γ ion 𝜌 𝑛 v 2 𝑛 ) + Γ rec 𝜌 𝑝 𝑇 𝑝 -Γ ion 𝜌 𝑛 𝑇 𝑛 𝛾 ( 𝛾 -1 ) . (15)  \nBoth fluids are subject to the ideal gas law. The factor of 2 in equation (12) is included to account for the electron pressure in the plasma fluid.  \nThe two-fluid collisional coupling is mediated by 𝛼 𝑐 ( 𝑇 𝑛 , 𝑇 𝑝 , 𝑣 𝐷 ) (Draine 1986; Zank et al. 2018), that is defined as:  \n𝛼 𝑐 = 𝛼 𝑐 ( 0 ) √︂ 𝑇 𝑛 + 𝑇 𝑝 2 √︄ 1 + 9 𝜋 64 𝛾 2 ( 𝑇 𝑛 + 𝑇 𝑝 ) 𝑣 2 𝐷 , (16)  \nwhere 𝛼 𝑐 ( 0 ) is the initial coupling and 𝑣 𝐷 = | v 𝑛 -v 𝑝 | is the magnitude of the drift velocity between the neutral components and the hydrogen plasma. The collisional frequencies are determined by the product of 𝛼 𝑐 and the fluids density: the ion-neutral collisional frequency is defined as 𝛼 𝑐 𝜌 𝑛 = 𝜏 -1 𝑖𝑛 , and the neutral-ion collisional frequency is 𝛼 𝑐 𝜌 𝑝 = 𝜏 -1 𝑛𝑖 , where 𝜏 𝑖𝑛 and 𝜏 𝑛𝑖 are the ion-neutral and neutral-ion coupling time scales, respectively. While the coupling between ions and neutrals is governed by 𝛼 𝑐 , the collisional coupling between ions and electrons is modelled by imposing a diffusivity 𝜂 in the system. The distribution for 𝜂 chosen in this work is spatially irregular and time-dependent: more details can be found in Section 2.1. Initially, the two fluids are in thermal and ionisation equilibrium.  \nThe terms Γ rec and Γ ion are the recombination and collisional ionisation rates for a hydrogen atom. The normalised empirical forms of the rates (Voronov 1997; Smirnov 2003) are:  \nΓ rec = 𝜌 𝑝 √︁ 𝑇 𝑝 √︁ 𝑇 𝑓 𝜉 𝑝 ( 0 ) 𝜏 IR , (17)  \nΓ ion = 𝜌 𝑝 𝑒 -𝜒 𝜒 0 . 39 0 . 232 + 𝜒 ! ˆ 𝑅 𝜉 𝑝 ( 0 ) 𝜏 IR , (18)  \n𝜒 = 13 . 6 𝑇 𝑓 𝑇 𝑒 0 𝑇 𝑝 , (19)  \nˆ 𝑅 = 2 . 91 · 10 -14 2 . 6 · 10 -19 √︁ 𝑇 𝑒 0 . (20)  \nTwo characteristic temperatures appear in Equations (17)-(20), and are based on a physical reference electron temperature 𝑇 0, expressed in Kelvin. 𝑇 𝑒 0 is the value of 𝑇 0 converted in electron volts. In our model we assume that protons and electrons are at the same temperature. Therefore, the choice of 𝑇 0 at the beginning of calculations determines the initial ion fraction 𝜉 𝑝 ( 0 ) . 𝑇 𝑓 is a normalisation factor defined as  \n𝑇 𝑓 = 𝛽𝛾 4 2 𝜉 𝑝 ( 0 ) 𝜉 𝑛 ( 0 ) + 2 𝜉 𝑝 ( 0 ) , (21)  \nand ensures that the reference temperature 𝑇 0 becomes the initial non-dimensional temperature of the simulation. The free parameter 𝜏 IR determines the ratio between recombination time scale and dynamic time: following our normalisation, 𝜏 IR states the initial recombination rate.  \nThe terms Φ 𝐼 and 𝐴 heat in equation (8) are associated to the ionisation potential and account for radiative losses. Φ 𝐼 approximates the energy removed by the system through ionisation, while 𝐴 heat is an arbitrary heating term included to obtain an initial equilibrium. Their non-dimensional forms are:  \nΦ 𝐼 = Γ ion 𝜌 𝑛 ˆ Φ , (22)  \n𝐴 heat = Γ ion ( 𝑡 = 0 ) 𝜌 𝑛 ( 𝑡 = 0 ) ˆ Φ , (23)  \nwhere ˆ Φ = 13 . 6 𝛽 2 𝐾 𝐵 𝑇 0 ensures consistency between the normalisation of ionisation potential and the equations modelling the system, which are normalised to the total Alfv\'en speed 𝑣 𝐴 = 1 and the total density 𝜌 0. As we include these heating/cooling terms associated to ionisation and recombination processes, the total energy of the system is not conserved in PIP. It is however conserved in fully-ionised MHD, as the solver is known to conserve energy in the absence of losses. More details on the atomic internal structure used to estimate the ionisation potential can be found in Snow & Hillier (2021).  \nThe physical variables in our system are non-dimensionalised as follows:  \nr → 𝑟 \' ˜ r , B → 𝐵 \' ˜ B , v → 𝑣 \' ˜ v , 𝑡 → 𝑡 \' ˜ 𝑡, 𝜌 → 𝜌 \' ˜ 𝜌, 𝑃 → 𝑃 \' ˜ 𝑃, (24)  \nwhere 𝑟 \' is the radius of the flux tube, 𝐵 \' = 1 is the initial axial field at 𝑟 = 0, 𝑣 \' = 𝑣 𝐴 = 𝐵 \' / √ 𝜌 \' as the permeability is set as 𝜇 0 = 1, the characteristic time scale 𝑡 \' = 𝑡 𝐴 = 𝑟 \' / 𝑣 𝐴 is the Alfv\'en transit time across the flux rope, 𝜌 \' = 1 is the total density and the total gas pressure 𝑃 \' = 𝛽 / 2. This non-dimensionalisation is consistent with the work in Hood et al. (2009).', '2.1. Initial conditions': "The initial setup is provided by a twisted magnetic flux tube in force-free equilibrium that is unstable to an ideal MHD kink instability. The initial non-potential, force-free magnetic field satisfies the condition:  \n∇ × B = 𝛼 ( r ) B , (25)  \nwhere 𝛼 = 𝜇 0 J · B / 𝐵 2 is the initial stress of the flux tube, and can relax through reconnection and reach a lower energy state, while preserving magnetic helicity. The field relaxes towards the minimum energy state, which must have a linear or constant 𝛼 -field profile: for this mechanism to work, the magnetic field must be sufficiently stressed beyond the relaxed state, before an instability releases the stored magnetic energy.  \nIn our simulations we model a cylindrical flux rope with a radius 𝑟 = 1. The magnetic field components, consistent with the case studied in Hood et al. (2009) for an initial smooth 𝛼 -profile, are reported in the following equations in cylindrical coordinates. For 𝑟 < 1 we set:  \n𝐵 𝜃 ( 𝑟 ) = 𝜆𝑟 ( 1 -𝑟 2 ) 3 , (26)  \n𝐵 𝑧 ( 𝑟 ) = √︂ 1 -𝜆 2 7 + 𝜆 2 7 ( 1 -𝑟 2 ) 7 -𝜆 2 𝑟 2 ( 1 -𝑟 2 ) 6 , (27)  \n𝛼 = 2 𝜆 ( 1 -𝑟 2 ) 2 ( 1 -4 𝑟 2 ) 𝐵 𝑧 ( 𝑟 ) , (28)  \nwhile for 𝑟 ≥ 1 these components become:  \n𝐵 𝜃 = 0 , 𝐵 𝑧 = √︂ 1 -𝜆 2 7 , 𝛼 = 0 . (29)  \nFigure 1. Left: initial 2D profile of 𝐵 𝑧 at 𝑧 = 0 (top panel) and 𝐽 𝑧 at 𝑧 = 0 (bottom panel) in the 𝑥𝑦 -plane for all simulations. The black contour lines in the bottom panel represent the field lines of the 𝐵 𝑥 𝐵 𝑦 vector field. Right: initial 3D setup of the twisted magnetic field lines of case M1. The slices at 𝑧 = -10 , 0 , 10 represent the magnitude of the plasma temperature.  \n<!-- image -->  \nThis choice of the initial magnetic field components leads to a configuration with zero net-current. In the equations above, 𝜆 is a constant parameter corresponding to a measure of the twist in the field. Due to the requirement from equation (27) that 𝐵 2 𝑧 must be positive, the value of 𝜆 is constrained to be 𝜆 < 64 / 965 √ 1351 = 2 . 438. In this work we impose 𝜆 = 1 . 8, consistent with Hood et al. (2009): the choice of these initial conditions also leads to a change of sign for 𝛼 at 𝑟 = 0 . 5. Figure 1 shows the initial profiles of 𝐵 𝑧 and the current density component 𝐽 𝑧 in the 𝑥𝑦 -plane ( 𝑧 = 0) set for all simulations. The current density profile is continuous across the domain, and a smooth transition occurs from 𝐽 ≠ 0 to 𝐽 = 0 at 𝑟 = 1. The three-dimensional form of the magnetic field can be observed in Figure 1, compared to the initial plasma temperature of case M1.  \nThe plasma 𝛽 associated with the bulk pressure is 0.1, consistent with the upper chromosphere (Gary 2001). Chromospheric neutrals play an important role in the decay of the non force-free components of the photospheric field, as it transitions to the force-free coronal field (Goodman 2000; Khodachenko et al. 2004; Arber et al. 2009). In a two-fluid system, the balance with the Lorentz force can only be achieved by varying the plasma pressure: as the charges have an effective 𝛽 of 0.02, as such we expect the system to be approximately force-free. Therefore, a force-free equilibrium is a suitable starting point for this work.  \nThe calculations are run with a non-uniform, current-dependent resistivity, that accounts for electron-ion collisions and is determined as follows:  \n𝜂 = 𝜂 0 +       0 , | 𝐽 | < 𝐽 crit 𝜂 1 , | 𝐽 | > 𝐽 crit (30)  \nwhere 𝜂 0 is the background uniform resistivity, set equal to 𝜂 0 = 10 -4 , and 𝜂 1 = 10 -3 is the anomalous resistivity component that is only switched on when the magnitude of the current exceeds a critical value, chosen to be 𝐽 crit = 5 to be above the initial maximum current magnitude determined by the initial conditions. This choice of diffusivity, resolved by the grid, is similar to the work of Arber et al. (1999), also used by Browning et al. (2008) and Hood et al. (2009). As ambipolar diffusion represents strongly coupled two-fluid effects in a single-fluid regime (Khomenko & Mart'ınez-G'omez 2024), its role is already covered by the more complex two-fluid physics of our setup. As such, no separate resistive term (i.e. ambipolar diffusion, Brandenburg & Zweibel 1994, 1995; Chae & Litvinenko 2021) is included to account for the effects of neutral-charge interactions.  \nIn order to trigger the instability in the system, a small velocity perturbation is included. The initial velocity perturbation, imposed in both the ion and neutral fluids to break the initial equilibrium, is a 𝑥𝑦 -plane white noise of the form:  \n𝑣 𝑥, 𝑝 , 𝑣 𝑦, 𝑝 = 0 . 05 𝑒 -( 𝑥 2 + 𝑦 2 ) · ( random noise ) . (31)  \nThe perturbation results in a larger velocity magnitude at the centre of the flux rope, perpendicular to the flux rope axis, and a smaller velocity magnitude at the extremities of the flux rope. As this perturbation is set to be the same for all cases, its effects  \nTable 1. List of the simulation parameters.  \nmight vary following the different coupling between ions and neutrals, resulting in the onset of different modes. Evidence of this process is discussed in Section 3.  \nThe simulations are run in a Cartesian computational domain 𝑥 = [-2 , 2 ] , 𝑦 = [-2 , 2 ] and 𝑧 = [-10 , 10 ] with a uniform grid. Although the analytical equilibrium is presented in cylindrical coordinates, the field components have been mapped onto the Cartesian grid to provide the initial states for the simulations. The domain is resolved by 496 × 496 × 868 grid points, corresponding to a cell size of Δ 𝑥 = 8 · 10 -3 , Δ 𝑦 = 8 · 10 -3 and Δ 𝑧 = 2 . 3 · 10 -2 . All boundaries are periodic, mimicking previous studies of instabilities in flux ropes (Shafranov 1957; Kruskal et al. 1958). Following this choice, the system simulates a potentially infinite number of infinitely long flux ropes, but the boundaries are set sufficiently far that the ropes never interact with each other, and the kink instability can be studied for a single structure.", '3. RESULTS': "We analyse the kink instability in three simulations, corresponding to a fully-ionised plasma case (M1) and two PIP cases (P1 and P2), run at different collisional coupling. M1 represents the limit for a completely coupled system ( 𝛼 𝑐 →∞ ), where neutrals and charges act as a single fluid. P1 corresponds to a weakly coupled ion-neutral plasma, where collisions occur on timescales ∼ 𝜌 -1 for a reference density 𝜌 , while P2 is run with a stronger collisional coupling between the fluids, with collisions occurring on timescales ∼ 0 . 1 𝜌 -1 . The parameter 𝜏 IR is varied in the PIP cases to maintain the same proportionality between collision frequencies and ionisation/recombination rates. The full array of parameters is shown in Table 1.  \nThe development of the kink instability is shown in Figure 2 for three evolution stages in both fully-ionised plasmas (M1) and partially-ionised plasmas (P1 and P2). At the onset of the instability, the current density is larger at the centre of the flux rope, as shown by panels ( 𝑎 ), ( 𝑏 ) and ( 𝑐 ). Under the distortion generated by the kink, the plasma flow pushes the field lines to form current sheets, that appear in panels ( 𝑑 ), ( 𝑒 ) and ( 𝑓 ) of Figure 2 as dark, thin and elongated regions where 𝐽 > 5. As the magnetic tension reduces with the straightening of the field lines, and the gas pressure increases with the thermal energy released from reconnection, the flux rope expands. During the expansion phase, shown at early times in panels ( 𝑑 ), ( 𝑒 ) and ( 𝑓 ), current sheets form at different locations and the reconnection process moves to regions further away from the centre of the flux rope. At later stages, shown in panels ( 𝑔 ), ( ℎ ) and ( 𝑖 ) reconnection occurs mostly near the external boundary of the flux rope for P1 and P2, while in case M1 several structures of high concentrations of current density are still present around the centre of the flux rope.  \nAs shown by Figure 2, where the same instability phase appears in the bottom panels at 𝑡 = 120 for M1, 𝑡 = 105 for P2 and 𝑡 = 65 for P1, respectively, the kink instability evolves much faster in partially-ionised plasmas than in fully-ionised plasmas. By comparing the two PIP cases, the flux rope expansion is faster when the collisional coupling 𝛼 𝑐 is smaller. The same development stage of the instability occurs earlier ( 𝑡 = 65) for case P1 (where 𝛼 𝑐 ( 0 ) = 1) than case P2 ( 𝑡 = 105, where 𝛼 𝑐 ( 0 ) is set to be equal to 10). This occurs as the growth of the kink instability scales with the Alfv'en speed, which is inversely proportional to the density. In case M1 the fluid is completely ionised, and the ion density coincides with the bulk density. In both PIP cases, run with the same bulk density as M1, the effective Alfv'en speed is larger as it depends on a combination of the smaller ion density and the density of the neutrals that are collisionally coupled to the charges: for this reason, the kink instability grows faster. The lower coupling in case P1 speeds up the reconnection rate and pushes the flux rope to expand faster than in case P2, where the stronger coupling between ions and neutrals results in the evolution of the system more closely resembling that of a single, completely ionised fluid.  \nThe charges temperature 𝑇 𝑝 and the current density magnitude at the three instability phases in Figure 2 are shown in Figure 3 in the 𝑥𝑧 -plane ( 𝑦 = 0) laying along the flux rope axis for cases P1 (top six panels), P2 (central six panels) and M1 (bottom six panels). The writhe of the flux rope axis is observed in the top two panels of each set, corresponding to 𝑡 P1 = 45, 𝑡 P2 = 75 and 𝑡 M1 = 90. In case P1, where the two fluids are almost decoupled, four crests located at regular distances from each other are observed along the 𝑧 -axis, which suggests the onset of a single, dominant mode of the instability. Similarly, four crests are  \nFigure 2. Contour plot of the current density magnitude 𝐽 at 𝑧 = 0, showing the development of the kink instability in M1 (left), P2 (centre) and P1 (right). Each row displays the same stage of development of the kink instability. Panels ( 𝑎 ), ( 𝑏 ) and ( 𝑐 ) show the onset of the kink instability. Panels ( 𝑑 ), ( 𝑒 ) and ( 𝑓 ) display the evolution at a later stage, and the formation on current sheets inside the flux rope. Panels ( 𝑔 ), ( ℎ ) and ( 𝑖 ) show the flux rope expanding and the formation of further current sheets. The arrows show the direction of the plasma flow v 𝑝 . Times are given in the same non-dimensional unit for all three cases. 𝐽 is saturated by fixing the maximum value of the colour scale to the critical value 𝐽 crit = 5, in order to enhance the structures characterised by strong currents.  \n<!-- image -->  \nobtained in case M1. In the intermediate case, however, it is possible to identify five crests, with a new, central peak appearing around 𝑧 = 0. While a dominant mode is observed both in the PIP case with the lowest coupling and in the full coupling limit (M1), the very presence of oscillations with crests located at irregular distances in P2 suggests the overlapping of two modes of the instability. Therefore, the partially coupled neutrals leads to the selection of different modes of the kink instability.  \nFigure 4 shows the time variation of the energy components integrated over the volume 𝑉 of the simulation box, and the maximum current density of cases M1, P1 and P2. In order to better show the variation of the energy terms, the net variation of internal and magnetic energy are displayed. The total kinetic energy is shown fully, and it differs from zero at 𝑡 = 0 due to the initial velocity perturbation.  \nThe passage from the linear growth phase to the nonlinear reconnection phase of the kink instability is identified by the decrease of magnetic energy Δ 𝑀 (panel 𝑏 of Figure 4), defined as:  \nΔ 𝑀 = ∫ 𝑉 𝐵 2 2 𝑑𝑉 -∫ 𝑉 𝐵 2 ( 𝑡 = 0 ) 2 𝑑𝑉, (32)  \nand the increase of kinetic and internal energy components (panels 𝑎 and 𝑐 of Figure 4) due to the energy conversion from reconnection, which for case M1 starts at 𝑡 ∼ 90, for case P1 is at 𝑡 ∼ 50 and for case P2 is at 𝑡 ∼ 80. The maximum current density reached in all three cases is approximately similar, as suggested by the trend of 𝐽 displayed in panel ( 𝑑 ).  \nt = 55  \nt = 105  \nt = 105  \nt = 65  \nt = 65  \nt = 65  \nt = 65  \nt = 90  \nt = 105  \nt = 105  \nt = 75  \nt = 75  \nt = 45  \nt = 45  \nt = 90  \nt = 90  \nt = 65  \nt = 105  \nt = 105  \nt = 90  \nt = 55  \nt = 45  \nt = 45  \nt = 75  \nt = 75  \nt = 90  \nt = 90  \nt = 105  \nt = 105  \nt = 105  \nt = 45  \nt = 55  \nt = 65  \nt = 55  \nt = 90  \nt = 75  \nt = 90  \nt = 90  \nt = 105  \nt = 120  \nt = 105  \nt = 120  \nt = 55  \nt = 65  \nt = 65  \nt = 90  \nt = 105  \nt = 105  \nt = 120  \nt = 90  \nt = 75  \nt = 105  \nt = 65  \nt = 120  \nt = 90  \nt = 75  \nt = 90  \nt = 75  \nt = 90  \nt = 105  \nt = 90  \nt = 75  \nt = 90  \nt = 105  \nt = 90  \nt = 105  \nt = 105  \nt = 120  \nt = 120  \nt = 90  \nt = 105  \nt = 105  \nt = 105  \nt = 120  \nt = 90  \nt = 105  \nt = 120  \nt = 90  \nt = 105  \nt = 90  \nt = 90  \nt = 105  \nt = 120  \nt = 105  \nt = 120  \nt = 105  \nt = 120  \nt = 120  \nt = 105  \nt = 105  \nt = 65  \nt = 65  \nt = 55  \nt = 55  \nt = 65  \nt = 65  \nt = 90  \nt = 90  \nt = 105  \nt = 105  \nt = 75  \nt = 75  \nt = 45  \nt = 45  \nt = 90  \nt = 90  \nt = 65  \nt = 65  \nt = 65  \nt = 105  \nt = 105  \nFigure 3. Contour plot of the plasma temperature 𝑇 𝑝 (left) and current density magnitude 𝐽 (right) in the 𝑥𝑧 -plane ( 𝑦 = 0), showing the development of the kink instability in case P1 (top panels), P2 (central panels) and M1 (bottom panels). Each set of panels displays the same stage of development of the kink instability shown in Figure 2. 𝑇 𝑝 is saturated by fixing the maximum value of the colour scale to the value 0.1, to enhance the plasma structures where the largest temperature increase occurs. Similarly, 𝐽 is saturated by fixing the maximum value of the colour scale to the critical value 𝐽 crit = 5 in order to enhance the structures characterised by strong currents.  \n<!-- image -->  \nRegarding the distribution of the converted magnetic energy, the individual contributions of plasma (blue) and neutrals (red) to the kinetic energy in panel ( 𝑎 ) and the total kinetic energy, defined as:  \n𝐾 TOT = ∫ 𝑉 1 2 𝜌 𝑝 v 2 𝑝 + 1 2 𝜌 𝑛 v 2 𝑛 ! 𝑑𝑉, (33)  \nare smaller than the MHD case, as shown by the dashed black line (P1) and the dotted black line (P2), while the variation of the total internal energy 𝐼 TOT , shown in panel ( 𝑐 ) and calculated as:  \nΔ 𝐼 TOT = ∫ 𝑉 𝑝 𝑝 + 𝑝 𝑛 𝛾 -1 -𝑝 𝑝 ( 𝑡 = 0 ) + 𝑝 𝑛 ( 𝑡 = 0 ) 𝛾 -1 ! 𝑑𝑉, (34)  \nis larger in the PIP cases than in the MHD case (solid black line) when compared at the same growth phase of the instability, despite smaller single contribution from charges and neutrals.  \nt = 55  \nt = 90  \nFigure 4. Kinetic energy ( 𝑎 ), time variation of ( 𝑏 ) magnetic energy Δ 𝑀 and ( 𝑐 ) internal energy Δ 𝐼 , and maximum current density ( 𝑑 ) for M1, P1 and P2. In all four panels, the solid lines are associated to variables of the MHD case (M1), the dashed lines refer to case P1 and the dotted lines refer to case P2. In panels ( 𝑎 ) and ( 𝑐 ) blue lines are associated to the plasma energy components, red lines are associated to the neutral energy components and black lines are associated to the bulk energy component (only plasma for case M1, plasma and neutrals for cases P1 and P2). The kinetic energy is the only energy component that is perturbed initially.  \n<!-- image -->  \nThe difference in the increase of the bulk internal energy between PIP and MHD cases means that the temperature of cases P1 and P2 grows faster. In general, it must be expected that the temperature is higher in the PIP cases than in the MHD case: this comes from the inclusion of further heating terms that result from the interaction of charges and neutrals. This will be discussed in more detail in Section 3.2.  \nFollowing the work of Browning et al. (2008), the instability growth rate can be calculated as:  \n𝜎 = 1 2 𝑑 log ( 𝐾𝐸 ) 𝑑𝑡 . (35)  \nTherefore, in this work, the linear growth rate of the kink instability is estimated as the angular coefficient of the slope given by the initial increase of the natural logarithm of the total kinetic energy (black lines in panel 𝑎 Figure 4), divided by a factor of 2. The dimensionless growth rate of the three cases is, respectively,  \n- · 𝜎 (P1) = 0 . 18 ± 0 . 02, in the interval 𝑡 = 20 -60,\n- · 𝜎 (P2) = 0 . 118 ± 0 . 006, in the interval 𝑡 = 30 -75,\n- · 𝜎 (M1) = 0 . 082 ± 0 . 004, in the interval 𝑡 = 55 -115,  \nwhere the error is determined as the standard deviation from the linear fit of the kinetic energy logarithm.  \nThe estimated growth rate decreases at the increase of 𝛼 𝑐 . The variation in growth rate between the fully coupled case (M1) and the least collisionally coupled case (P1) is of a factor of 2.2. The larger growth rate comes as a consequence of the variation in effective plasma density: in case of a fully decoupled case ( 𝛼 = 0) run at the PIP ion fraction = 0 . 1, the growth of the instability is expected to be a factor ∼ √ 10 (∼ 3 . 16 ) larger than the MHD case presented here. The value of the growth rate for our intermediate coupling cases prove that there is a direct response of the coupling in affecting the instability growth rate.", '3.1. Two-fluid effects': "In a partially-ionised plasma, the interaction of ions and neutrals can lead to further dynamics that might affect the physics of the kink instability. Here case P1 is examined at 𝑡 = 90, a late stage of the nonlinear evolution of the kink instability. Figure 5 shows the temperature difference 𝑇 𝑛 -𝑇 𝑝 , the magnitude of the drift velocity between fluids 𝑣 𝐷 = | v 𝑛 -v 𝑝 | , the ionisation rate  \nFigure 5. Contour plot of ( 𝑎 ) the temperature difference 𝑇 𝑛 -𝑇 𝑝 , ( 𝑏 ) the drift velocity magnitude | 𝑣 𝐷 | , ( 𝑐 ) the current density magnitude 𝐽 , ( 𝑑 ) the ionisation rate Γ ion and ( 𝑒 ) the recombination rate Γ rec for the PIP case P1 at 𝑡 = 90 (top panels), and details of the same variables in the area around a current sheet (bottom panels). The variables are shown at the centre of the flux rope ( 𝑧 = 0). The magnitude of ionisation and recombination rates is presented with a logarithmic color scale. The area covered by the bottom panels corresponds to the region identified in the black boxes in the top panels.  \n<!-- image -->  \nΓ ion, the recombination rate Γ rec and the current density magnitude 𝐽 at the centre of the flux rope ( 𝑧 = 0). The whole section of the flux rope is shown in the top panels, while the detail of a smaller area surrounding a current sheet is displayed in the bottom panels.  \nThe temperature distribution differs between neutral and plasma, as shown by the panels in column ( 𝑎 ) of Figure 5 where golden areas represent the regions where the plasma is hotter than the neutrals and purple areas identify the regions where the neutral fluid is hotter than the plasma. The plasma temperature exceeds the neutral temperature mostly in regions of high current density magnitude, as shown by the comparison between the panels in column ( 𝑎 ) and ( 𝑐 ), while patches of higher neutral temperature than the plasma are observed outside the reconnecting regions. The maximum temperature difference reached in these areas ( ∼ 0 . 17) is very large when compared to the background temperature of 0.075 of both species, and dimensionally is equivalent to a difference of about 3 · 10 4 K. The magnitude of the drift velocity is also very high, reaching values of about a fifth of the bulk Alfv'en speed, mostly in correspondence of high current density structures as shown by the comparison with 𝐽 . This means that in high current density regions the two fluids are mostly decoupled, a feature that is also shown by the large difference in temperature in these areas, with peaks in plasma temperature within the current sheets and peaks in neutral temperature immediately outside.  \nIn P1, ionisation and recombination processes occur with a comparable increase inside the flux rope. Initially Γ ion = 10 -6 and Γ rec = 10 -5 , values that remain constant in the background plasma. Spikes of ionisation rate are observed within the areas of larger current density and 𝑇 𝑝 , with a maximum value Γ ion = 2 . 4 · 10 -3 at 𝑡 = 90: inside the current sheet, the plasma compression and high temperature allow to ionise the neutrals dragged by collisions with the ions. The larger ionisation rate leads to a thickening of the current sheet, as noted by Murtas et al. (2022). Outside the longer current sheets, recombination rates increase proportionally to the plasma temperature, as shown by the comparison with the temperature difference map in column ( 𝑎 ), with peaks of Γ rec = 1 . 2 · 10 -5 at 𝑡 = 90. A similar trend is found for the case with a larger ion-neutral coupling (P2). In case P2, the background Γ ion = 10 -5 and Γ rec = 10 -4 , with peaks of Γ ion = 4 . 4 · 10 -2 and Γ rec = 1 . 4 · 10 -4 at 𝑡 = 105.  \nAs shown by the zoom-in of the area around a current sheet, displayed in the bottom panels of Figure 5, peaks of drift velocity (column b ) occur in the outflow regions of the current sheet, while in the inflow both neutral and plasma velocity are comparable. This has previously been observed in the results presented in Murtas et al. (2021) and Murtas et al. (2022), where an evident  \ndecoupling of the species occurred in the outflow regions of current sheets. The expulsion of neutrals from the current sheet lead to a more turbulent mixing in the regions around, and to local peaks of the neutral temperature.  \nFigure 6 shows the time evolution of the kink instability and the variation of the drift velocity at three locations of the flux rope. At all times the drift velocity produces similar structures of comparable magnitude along the flux rope. In the early stages of the simulation ( 𝑡 = 0) characterised by the largest writhe, the drift velocity is negligible. The first structures in 𝑣 𝐷 appear when the nonlinear phase of the instability begins ( 𝑡 = 30), and their distribution along a spiral structures follows the relaxation of the magnetic field lines. At 𝑡 = 60 the drift velocity increases are concentrated primarily at the external surface of the flux rope, where the first current sheets are formed, as expected from the distribution of 𝑣 𝐷 magnitude shown in Figure 5 and discussed above. Finally, the velocity structures created by the drift increase in complexity with the increasing turbulent motions occurring within the flux rope at later times ( 𝑡 = 90).", '3.2. Temperature and heating terms': "Beyond the term of Ohmic heating, the high velocity drifts observed in partially-ionised plasmas and discussed in Section 3.1 can lead to heating generated by collisions taking place between ions and neutrals at these high speeds. In this Section we examine the evolution of the fluids temperature and of the heating terms with the growth of the kink instability.  \nThe time variation of both neutral and plasma temperatures are shown in Figure 7 for case P1 and Figure 8 for case P2 at the same evolution phases displayed in Figure 2. In a partially-ionised plasma, fluids can be heated through Ohmic heating 𝜂𝐽 2 , generated by the current dissipation in regions at high current densities, and by frictional heating, which is generated by the drift between ions and neutrals. Frictional heating is characterised by two different components (Snow et al. 2023). The first component is a collisional term related to the drift velocity between ions and neutrals. The non-dimensional form of the collisional frictional heating is:  \n𝐹 heat , 1 = 1 2 𝛼 𝑐 ( 𝑇 𝑛 , 𝑇 𝑝 , 𝑣 𝐷 ) 𝜌 𝑛 𝜌 𝑝 𝑣 2 𝐷 , (36)  \nand depends on both the collisional frequencies and the drift velocity between charges and neutrals. The second term adding to frictional heating is linked to the work done through ionisation-recombination processes, and takes the form:  \n𝐹 heat , 2 = Γ rec 𝜌 𝑝 v 2 𝑝 - ( Γ rec 𝜌 𝑝 + Γ ion 𝜌 𝑛 ) v 𝑛 · v 𝑝 + Γ ion 𝜌 𝑛 v 2 𝑛 . (37)  \nThe contour maps of the Ohmic heating and the collisional frictional heating terms are shown in Figure 7 for case P1 and 8 for case P2. In both Figures, the colour scales are set equal for an easier comparison between the two PIP cases. At the beginning of the kink instability both 𝑇 𝑝 and 𝑇 𝑛 of cases P1 and P2 present very small variations with respect to the initial conditions, as shown by the panels in row ( 𝑎 ) of Figure 7-8. The plasma temperature starts increasing when the first current sheets form inside the flux rope, as shown in rows ( 𝑏 ) and ( 𝑐 ) of Figure 7-8, and as a consequence of the thermal coupling between the two fluids the neutral temperature starts to increase in the same areas of high 𝑇 𝑝 . The Ohmic term is larger in magnitude than frictional heating, as shown by the comparison of the heating terms magnitude in both Figure 7-8, and its action is localised to the current sheets structures. Outside these areas, collisional frictional heating provides a comparable contribution to heating both charges and neutrals. The frictional heating produced by ionisation-recombination processes is a few orders of magnitude smaller in both PIP cases, reaching a peak of 1 . 78 · 10 -5 for P1 at 𝑡 = 65 and 4 . 13 · 10 -5 for P2 at 𝑡 = 105, hence its contribution can be considered negligible.  \nThe temperatures and heating terms of the two PIP cases can be compared at similar growth stages of the instability, and across the central section of the flux rope. In case P1, where the initial collisional coupling is 𝛼 𝑐 ( 0 ) = 1, the plasma temperature at the later stages displays the largest increase, reaching peaks of 𝑇 𝑝,𝑃 1 ∼ 0 . 25, 3.3 times the background temperature in the reconnecting regions (corresponding to peaks of ∼ 4 . 2 · 10 4 K), against 𝑇 𝑝,𝑃 2 ∼ 0 . 21 for case P2 where 𝛼 𝑐 ( 0 ) = 10, which is about 2.7 times the background temperature (corresponding to ∼ 3 . 4 · 10 4 K). From the underlying empirical ionisation and recombination equations 17-18, the equilibrium ionisation fraction can be defined as a function of temperature only (Snow & Hillier 2021), i.e.,  \n𝜉 𝑝 = 1 𝐹 ( 𝑇 ) 𝐺 ( 𝑇 ) + 1 , (38)  \nwhere 𝐹 ( 𝑇 ) , 𝐺 ( 𝑇 ) are the temperature dependent components of the recombination and ionisation rates respectively. The temperature maxima in both P1 and P2 cases is high enough that the plasma should be, locally, completely ionised ( 𝜉 𝑝 ≈ 1) based on Equation 38. However, Equation 38 is based on steady-state equilibrium, whereas the simulation is time-dependent and allows departures from equilibrium to exist. The ionisation and recombination rates are relatively slow, with Γ ion ≈ 10 -3 in  \nFigure 6. Time series of the kink instability of case P1 showing magnetic field lines (grey) and slices of the drift velocity squared ( v 𝑛 -v 𝑝 ) 2 . The times showed are, from top to bottom, 𝑡 = 0, 𝑡 = 30, 𝑡 = 60 and 𝑡 = 90.  \n<!-- image -->  \nFigure 7. Left to right: neutral temperature 𝑇 𝑛, norm and plasma temperature 𝑇 𝑝, norm normalized by the respective background temperatures, collisional frictional heating 𝐹 heat , 1 and Ohmic heating 𝜂𝐽 2 of case P1 at 𝑡 = 45 (row 𝑎 ), 𝑡 = 55 (row 𝑏 ) and 𝑡 = 65 (row 𝑐 ).  \n<!-- image -->  \nFigure 8. Left to right: neutral temperature 𝑇 𝑛, norm and plasma temperature 𝑇 𝑝, norm normalized by the respective background temperatures, collisional frictional heating 𝐹 heat , 1 and Ohmic heating 𝜂𝐽 2 of case P2 at 𝑡 = 75 (row 𝑎 ), 𝑡 = 90 (row 𝑏 ) and 𝑡 = 105 (row 𝑐 ).  \n<!-- image -->  \nthe reconnection region, as shown in Figure 5 for case P1 where Γ ion peaks at 2 . 4 · 10 -3 at 𝑡 = 90. As such, the timescale for ionisation-recombination equilibrium to be obtained is long compared to the life of the feature. Within the reconnection region, the ionisation fraction is roughly 0.1 and as such, the local neutral density is far larger than would be predicted by ionisationrecombination equilibrium. The larger temperature increase in P1 (the case with the smallest ion-neutral coupling) depends on the difference in the effective plasma 𝛽 , which is smaller in the most decoupled case as a lower coupling results in a lower gas density. The largest increases in neutral temperatures for the two PIP cases are closer, being 𝑇 𝑛,𝑃 1 ∼ 0 . 11 and 𝑇 𝑛,𝑃 2 ∼ 0 . 13, 1.5 and 1.74 times the background temperature respectively, despite the difference in the magnitude of the plasma temperature.  \nFigure 9. Global Ohmic heating (black) and collisional frictional heating (red) as a function of time for case M1 (solid line), P1 (dashed line) and P2 (dotted line).  \n<!-- image -->  \nAsecond difference is observed in the magnitude of the heating components between case P1 and P2. While the Ohmic heating 𝜂𝐽 2 is similar in both cases, the maximum being ∼ 0 . 25 for case P1 and ∼ 0 . 26 for case P2 in the last phase of the instability (row 𝑐 of Figure 7-8), the contribution of the collisional frictional heating is larger in case P2, being ∼ 0 . 04 against a maximum value of ∼ 0 . 014 for case P1. The difference in frictional heating comes from the difference in the collision coupling between the two cases. Collisional frictional heating is directly proportional to 𝛼 𝑐 as shown by the definition in equation (36), and directly affects the increase in temperature of both charges and neutrals.  \nWhile the collisional frictional heating is smaller locally than the Ohmic heating, its total contribution is larger than Ohmic heating in both PIP cases. This is shown in Figure 9, where Ohmic heating and collisional frictional heating integrated over the domain are displayed as a function of time. This occurs because the Ohmic heating, despite being larger in magnitude, is localised in the few current sheet structures at higher current densities, while the frictional heating produced by the drift between charges and neutrals is produced more uniformly inside the whole flux rope. As already noticed from the magnitude of the collisional frictional heating in Figure 7-8, P2 develops larger amounts of collisional frictional heating than P1 due to the stronger coupling between fluids.  \nFigure 10 shows the changes of the mean plasma and neutral temperatures, averaged over the whole domain, of all three simulations. Due to the normalization based on the Alfv'en speed, the initial mean plasma temperature is slightly higher in the MHD case, where 𝑇 𝑝 ( 𝑡 = 0 ) = 0 . 0833, than in the PIP cases, where at 𝑡 = 0 𝑇 𝑝 = 𝑇 𝑛 = 0 . 0758. The time variation of the mean temperature is very small, due to averaging across a large portion of the domain where the plasma is unperturbed. Comparing similar stages of the instability, the plasma temperature of M1 increases up to 𝑇 𝑝 = 0 . 0842 with a rise of 1% at 𝑡 = 120, and a similar increase is seen for P2, where 𝑇 𝑝 = 0 . 0768 (1.3%) at 𝑡 = 105, while for case P1 𝑇 𝑝 = 0 . 077 (1.6%) at 𝑡 = 65. Due to the original difference in background temperature, the mean temperature of case M1 is higher than P1 and P2, but the PIP cases show a temperature increase that is inversely proportional to the two-fluid coupling.", '4. SUMMARY AND DISCUSSION': "The kink instability of flux ropes may take place in several layers of the solar atmosphere, but has been well modelled only in coronal plasmas. Due to the lack of numerical models, it is still unknown of how this instability develops in a plasma that is only partially-ionised, such as the chromospheric plasma. In this work we present a study of the kink mode developing in PIP, as the growth of the instability is observed in presence of partial ionisation and compared to the evolution in a fully-ionised plasma. Our results can be summarised as follows.  \n- · The growth time scale of the helical kink instability increases with the collisional coupling. The nonlinear phase begins earlier in PIP, and is faster for the most decoupled case (P1). This depends on the scaling with respect to the Alfv'en speed, which is larger the lower the ion fraction and the coupling with neutrals. Consequently, the growth rate is larger the lower the ion-neutral collisional coupling. In all three cases, the estimated growth rate of the instability, normalised by the Alfv'en time, matches the values found in Browning et al. (2008) ( 𝜎 ∼ 0 . 05 -0 . 20) and Hood et al. (2009) ( 𝜎 ∼ 0 . 075). In the limit of a completely decoupled system ( 𝛼 𝑐 = 0), as pointed out by Murtas et al. (2021, 2022), we expect the instability of  \nFigure 10. Evolution of the global  \n<!-- image -->  \nmean plasma temperature (blue) and mean neutral temperature (red) of cases M1 (solid line), P1 (dashed line) and P2 (dotted line).  \nthe growth rate to be the fastest: this can be considered as a fully-ionised plasma, whose total density is equal to the density of the partially-ionised plasma.  \n- · The loss of magnetic energy during reconnection results in a larger increase of kinetic energy for the MHD case and of internal energy for the PIP cases. One of the reasons for the faster temperature increase of the PIP cases is that there are more heating terms, such as frictional heating, contributing.\n- · The maximum Ohmic heating is similar in both PIP cases and locally bigger than frictional heating, but frictional heating is more distributed and its total value in the system is larger. In PIP cases, collisional frictional heating is the major contribution to heating. In the weakly coupled regimes explored in this work, collisional frictional heating is bigger in case P2 than P1, exhibiting an increase that is directly proportional to the collisional coupling. In more strongly coupled systems, collisional frictional heating is expected to reduce, following the increasingly smaller drift velocity between charges and neutral.  \nThese results suggest that flux ropes forming in partially-ionised plasmas are more unstable to the kink mode, compared to fully-ionised plasmas of equal total density. This could lead to a faster release of energy following the shorter time scale of the instability, and contribute to the onset of chromospheric jets, where plasmoid formation has been already observed as a product of reconnection (Singh et al. 2011) and chromospheric micro- and mini-filament eruptions (Sterling et al. 2015, 2016; Sterling & Moore 2016; Samanta et al. 2019; Sterling et al. 2020). Using the values of temperature and plasma 𝛽 imposed as initial conditions, the Alfv'en speed is ∼ 40 km s -1 . The typical length of a mini-filament is ∼ 8 · 10 3 km (Sterling et al. 2015), therefore we estimate the initial flux rope radius to be ∼ 400 km. This leads to a time scale 𝜏 𝐴 ∼ 10 s. For typical chromospheric parameters, including a total number density 𝑛 = 10 20 m -3 , the ion-neutral collisional frequency is calculated as  \n𝜈 𝑖𝑛 = 𝑛 𝑛 √︂ 8 𝐾 𝐵 𝑇 𝜋𝑚 𝑖𝑛 𝜎 𝑖𝑛 , (39)  \nwhere 𝑛 𝑛 is the neutral number density is a fraction of the total density, 𝑚 𝑖𝑛 = 𝑚 𝑖 𝑚 𝑛 /( 𝑚 𝑖 + 𝑚 𝑛 ) ∼ 𝑚 𝑖 and 𝜎 𝑖𝑛 = 5 · 10 -19 m 2 is the cross section, and varies in the interval 10 3 -10 6 s -1 (Leake et al. 2005; Khomenko & Collados 2012). In non-dimensional form, this would correspond to an 𝛼 𝑐 ∼ 10 4 -10 7 , which is a strongly coupled regime. The effects of weak ionisation and coupling are enhanced in much smaller structures: for flux ropes with radius on the order of 10s of meters the dimensional ion-neutral collisional frequency 𝜈 𝑖𝑛 = 0 . 09 -0 . 9 s -1 . This however has to be considered as a lower limit for which strong two-fluid effects can be observed: both ion fraction and plasma 𝛽 can be orders of magnitude lower than the values imposed here, which would lead to stronger local two-fluid effects. It has to be mentioned that in this work we haven't explored the regime between weakly and strongly coupled plasmas ( 𝛼 𝑐 = 10 2 -10 3 ), which would lead to the enhancement of partial ionisation effects in larger structures. We defer this to future studies.  \nPartial ionisation is also responsible for more heating than what produced during the instability growth in fully-ionised plasmas. The new heating term is generated by the drifts between charges and neutrals (frictional heating) that is dominant over the Ohmic heating. Flux ropes are therefore more heated in partially-ionised plasmas, as more energy is converted into heat than kinetic energy, unlike the fully-ionised coronal processes. This result identifies an important contribution to the heating of the chromosphere, and can be connected to the works of Gudiksen & Nordlund (2005) and Hansteen et al. (2007), whose simulations suggest that most of the heating in the solar atmosphere occurs at chromospheric heights. The localised heating, produced by the small-scales explosive events that are triggered by the instability of chromospheric flux ropes (e.g. chromospheric jets), provides an important contribution to the chromospheric heating, and could also significantly fuel the coronal heating (De Pontieu et al. 2009).  \nThe study of mini-filaments of finite length would require a different type of boundary conditions at the footpoints (i.e. conductive walls, Hood et al. 2009, line-tying conditions, Hood & Priest 1979). While periodic boundaries are still appropriate to examine how partial ionisation changes the evolution from the standard MHD case. This allows us to infer that even with more realistic boundaries we will see faster growth of the instability and more heating. We defer the use of a more realistic setup to future studies.", 'ACKNOWLEDGEMENTS': "The authors thank Prof. Rony Keppens and Dr. Claire Foullon for the precious comments and suggestions that led to the completion of this study. GM is supported by the NASA Grant No. 80HQTR21T0087. AH and BS are supported by STFC Research Grant No. ST/R000891/1 and ST/V000659/1. For the purpose of open access, the authors have applied a 'Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising. This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure. This research also used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231.", 'DATA AVAILABILITY': 'The data that support the findings of this study are available from the corresponding author upon reasonable request. The (PIP) code is available at the following url: https://github.com/AstroSnow/PIP. Details of the code and equations are available in Hillier et al. (2016).', 'REFERENCES': 'Arber, T. D., Botha, G. J. J., & Brady, C. S. 2009, ApJ, 705, 1183, doi: 10.1088/0004-637X/705/2/1183  \nArber, T. D., Longbottom, A. W., & Van der Linden, R. A. M.  \n1999, ApJ, 517, 990, doi: 10.1086/307222  \nBareford, M. R., & Hood, A. W. 2015, Philosophical Transactions of the Royal Society of London Series A, 373, 20140266, doi: 10.1098/rsta.2014.0266  \nBiskamp, D. 2000, Magnetic Reconnection in Plasmas  \nBotha, G. J. J., Arber, T. D., & Hood, A. W. 2011, A&A, 525, A96,  \ndoi: 10.1051/0004-6361/201015534  \nBrandenburg, A., & Zweibel, E. G. 1994, The Astrophysical  \nJournal Letters, 427, L91, doi: 10.1086/187372  \n-. 1995, The Astrophysical Journal, 448, 734,  \ndoi: 10.1086/176001  \nBrowning, P. K., Gerrard, C., Hood, A. W., Kevis, R., & van der  \nLinden, R. A. M. 2008, A&A, 485, 837,  \ndoi: 10.1051/0004-6361:20079192  \nChae, J., & Litvinenko, Y. E. 2021, Research in Astronomy and Astrophysics, 21, 232, doi: 10.1088/1674-4527/21/9/232  \nDe Pontieu, B., Carlsson, M., Rouppe van der Voort, L. H. M., et al. 2012, ApJL, 752, L12, doi: 10.1088/2041-8205/752/1/L12  \nDe Pontieu, B., McIntosh, S. W., Hansteen, V. H., & Schrijver, C. J. 2009, ApJL, 701, L1, doi: 10.1088/0004-637X/701/1/L1  \nDraine, B. T. 1986, MNRAS, 220, 133,  \ndoi: 10.1093/mnras/220.1.133  \nFan, Y. 2005, ApJ, 630, 543, doi: 10.1086/431733  \nFang, F., Fan, Y., & McIntosh, S. W. 2014, ApJL, 789, L19,  \ndoi: 10.1088/2041-8205/789/1/L19  \nGary, G. A. 2001, SoPh, 203, 71, doi: 10.1023/A:1012722021820  \nGoedbloed, J. P. H., & Poedts, S. 2004, Principles of  \nMagnetohydrodynamics  \nGoodman, M. L. 2000, ApJ, 533, 501, doi: 10.1086/308635  \nGordovskyy, M., Kontar, E. P., & Browning, P. K. 2016, A&A,  \n589, A104, doi: 10.1051/0004-6361/201527249'}