© ESO, 2016

1. Introduction

Solar flares are energetic phenomena commonly understood as fast releases of magnetic energy stored in the solar corona. The coronal plasma heated during such events produces characteristic emission signatures in the soft X-ray range, and the particles accelerated in the flaring regions produce non-thermal emission typically in the hard X-ray energy range (see e.g. the review by Benz 2008; Fletcher et al. 2011). Twisted magnetic flux ropes can store significant amounts of magnetic free energy in the solar corona and are susceptible to developing magneto-hydrodynamical (MHD) instabilities, which leads to the release of this energy. Coronal loops undergoing a kink instability, in particular, go through an initial exponential growth phase (the linear phase of the instability), until they start reconnecting with the background field (Browning et al. 2008). They then relax into a state with lower magnetic energy and less twist (hence releasing a fraction of the magnetic free energy that was initially stored). It has recently been shown that magnetic reconnection in helical (twisted) magnetic structures may result in impulsive plasma heating that is capable of producing the main properties of thermal X-ray emission in solar flares (Pinto et al. 2015) and in particle acceleration that is distributed within large volumes of plasma (including the lower corona and the chromosphere), which all leads to non-thermal X-ray emission (Gordovskyy et al. 2013, 2014). The kink instability may therefore be a viable alternative to the standard flare model, which assumes that particles accelerated to high energies in small volumes of plasma near the flaring loop apexes are transported downwards along the loop legs and then interact with the denser chromosphere to produce hard X-ray emission there as well as provide hot plasma that would thereafter fill up the loops and emit thermally in soft X-rays. The kink-instability scenario faces one long-standing difficulty: its triggering requires that the pre-flare loops are strongly twisted (see e.g. the discussion in Bareford et al. 2013), which observations of flaring coronal loops rarely seem to indicate (for both confined and ejective flares). Notable exceptions exist, however, such as reported by Srivastava et al. (2010), who measured a total twist angle of about 12π on a twisted coronal flux rope before the occurrence of a B5.0 class flare (see also simulations of this event by Botha et al. 2012). Pinto et al. (2015) pointed out, however, that the morphology of the thermal emission during the initial phases of a flare may lead to a strong underestimation of the actual twist of the underlying magnetic field lines. Furthermore, the kink-instability scenario, along with other scenarios resulting in distributed acceleration or re-acceleration (Turkmani et al. 2006; Cargill et al. 2006), may allow circumventing the so-called number problem (Brown 1976; Brown et al. 2009), which is one of the main difficulties of the standard flare model.

In the present study, we focus on deducing observational signatures inherent to the triggering of the kink instability in twisted coronal loops by combining thermal and non-thermal effects (plasma heating and particle acceleration) with application to confined flares. We consider here loop models consisting of twisted magnetic flux ropes embedded in a gravitationally stratified background (comprising a hot corona, a transition region, and a fraction of the chromosphere), which represent small (or moderate) confined flares. The origin of such twisted coronal loops can be thought of as the result of photospheric source rotation or (see e.g. Brown et al. 2003; Dalmasse et al. 2014) or of flux rope emergence through the photosphere (e.g. Luoni et al. 2011; Archontis et al. 2013; Pinto & Brun 2013). But the detailed analysis of the formation of such twisted coronal loops is beyond the scope of this study; we only discuss here the evolution of flux ropes that are already twisted and close to the onset of the kink instability.

We compare light curves and time-dependent intensity maps of thermal (free-free continuum) and non-thermal (bremsstrahlung caused by energetic electrons) radiation in the keV and deca-keV ranges. We investigate how the spatial distribution of the emitted photon fluxes relates to the dynamical and geometrical properties of the simulated loops, how the total continuum emission spectra evolves in time, and on how the properties of the emission measures respond to the plasma heating processes occurring in the magnetic loops.

2. Model description

We study the evolution of a twisted coronal loop embedded in a magnetised and gravitationally stratified background atmosphere. The coronal loops are twisted magnetic flux ropes with large-scale curvature (the loops are nearly semi-circular), and both footpoints are anchored in the chromosphere. The thermo-dynamical evolution of the system is calculated using fully compressible resistive MHD by means of three-dimensional numerical simulations. The MHD computations provide the basis for estimating the thermal emission by the plasma (in the soft X-ray range) and of the hard X-ray emission produced by charged test particles (electrons). The methods employed are described in detail in the following sections (Sects. 2.1, 2.2, and 2.3).

2.1. MHD equations and model parameters

The evolution of the twisted loops is determined by solving the Eqs. (1)–(4) in a 3D (x,y,z) domain taking into account atmospheric stratification, thermal conduction, and radiative losses.

The system is described in terms of magnetic field B, plasma velocity v, density ρ, and specific energy ϵ: $\begin{array}{ccc}\frac{\mathit{\partial \rho }}{\mathit{\partial t}}& \mathrm{=}& \\ \frac{\mathit{\partial }{v}}{\mathit{\partial t}}& \mathrm{=}& \\ \frac{\mathit{\partial ϵ}}{\mathit{\partial t}}& \mathrm{=}& \mathrm{-}\mathrm{(}{v}\mathrm{·}{\nabla }\mathrm{)}\mathit{ϵ}\mathrm{-}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}\mathit{ϵ}{\nabla }\mathrm{·}{v}\\ & & \\ \frac{\mathit{\partial }{B}}{\mathit{\partial t}}& \mathrm{=}& \end{array}$where η is plasma resistivity and the current density is j = ∇ × B/μ₀, and g is the gravitational acceleration, which is constant throughout the whole domain. Plasma pressure and temperature can be expressed as p = (γ−1)ρϵ and . Here γ is the ratio of specific heats of the plasma, is the average particle mass (which is half the proton mas m_p, corresponding to fully ionised hydrogen plasma). The electric field that is required for the test-particle modelling can be expressed as ${E}\mathrm{=}\mathrm{-}{v}\mathrm{\times }{B}\mathrm{+}\mathit{\eta }{j}\mathit{.}$(5)The S_cond, Q_joule, Q_visc, and L_rad terms in Eq. (3) account for variations of the specific energy due to thermal conduction, Joule heating, viscous heating, and radiative losses, respectively. Thermal conduction is defined as ${\mathit{S}}_{\mathrm{cond}}\mathrm{=}{\nabla }\mathrm{·}\mathrm{[}\mathrm{(}\mathrm{-}{\mathit{\kappa }}_{\mathrm{0}}{\mathit{T}}^{\mathrm{5}\mathit{/}\mathrm{2}}{B̂}\mathrm{·}{\nabla }{\mathit{T}}^{\mathrm{)}}{{B̂}}^{\mathrm{]}}\mathrm{+}{\nabla }\mathrm{·}\left[\mathrm{-}{\mathit{\kappa }}_{\mathrm{0}}{\mathit{T}}^{\mathrm{5}\mathit{/}\mathrm{2}}\frac{{\mathit{B}}_{\min }^{\mathrm{2}}}{{\mathit{B}}_{\min }^{\mathrm{2}}\mathrm{+}{\mathit{B}}^{\mathrm{2}}}{\nabla }\mathit{T}\right]\mathit{,}$(6)where B_min is some very small parameter. Hence, with a non-zero magnetic field (B ≫ B_min), we obtain Braginskii conductivity (κ₀T^5/2) along the magnetic field, which is negligible in the perpendicular direction, however, while with a negligible magnetic field (B ≪ B_min) it is κ₀T^5/2 and isotropic. The following characteristic dimensional factors are used: B₀ = 4.47 × 10^-3 T, ρ₀ = 2 × 10^-12 kg m^-3 and L₀ = 10⁶ m, yielding the characteristic values for Alfvén velocity v₀ = 2.82 × 10⁶ m s^-1, time t₀ = 0.36 s, temperature T₀ = 3.23 × 10⁸ K, and electric field E₀ = 1.26 × 10⁴ V m^-1. The Joule heating is Q_joule = η/ρj², where η is the local resistivity η = η_back + η_crit. The background resistivity η_back is constant throughout the simulation domain and represents classical resistivity. Its value is ${\mathrm{10}}^{-6}\sqrt{{\mathit{\mu }}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}}{\mathit{L}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{B}}_{\mathrm{0}}^{-1}$ (corresponding to the Lundquist number 10⁶), which means that this resistivity component does not make any noticeable contribution to the system evolution. The critical resistivity is defined as ${\mathit{\eta }}_{\mathrm{crit}}\mathrm{=}\{\begin{array}{c}\\ {\mathrm{10}}^{-3}\sqrt{{\mathit{\mu }}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}}{\mathit{L}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{B}}_{\mathrm{0}}^{-1}\mathit{,}& \mathit{j}\mathrm{\ge }{\mathit{j}}_{\mathrm{crit}}\\ \mathrm{0}\mathit{,}& \mathit{j}\mathit{<}{\mathit{j}}_{\mathrm{crit}}\mathit{,}\end{array}$(7)where ${\mathit{j}}_{\mathrm{crit}}\mathrm{=}{\mathit{j}}_{\mathrm{1}}\mathit{\rho }\mathit{/}{\mathit{\rho }}_{\mathrm{0}}\sqrt{\mathit{T}\mathit{/}{\mathit{T}}_{\mathrm{0}}}$. This component represents anomalous resistivity with the threshold similar to that for ion-acoustic instability, when the current drift velocity exceeds local thermal velocity. The value of j₁ constant is chosen so that before the kink instability η_crit = 0 everywhere in the domain. The chosen resistivity values have no exact physical justification, as in most other MHD models. However, the chosen values appear to be realistic, as they provide realistic reconnection times (tens of seconds). The radiative losses are calculated using the function from Klimchuk et al. (2008). The shock viscosity and corresponding heating rate used in our simulations are described in Arber et al. (2001). The adopted viscosity values are not representative of real viscosity in the coronal plasma, and the corresponding viscous heating rate is usually lower than the ohmic heating rate due to anomalous resistivity. However, high velocity gradients during the fast reconnection phase can sometimes result in substantial viscous heating, comparable to or even higher than ohmic heating. Ultimately, the viscous effects result in conversion of magnetic energy to internal energy of the plasma and can be considered as an additional resistivity effective at locations with high velocity gradients (see Bareford et al. 2016, for a more detailed description of this effect). The above Eqs. (1)−(4) along with Eqs. (5)−(7) are solved using the Lagrangian remap code LARE3D (Arber et al. 2001). We consider four models: three with small coronal loops with different magnetic field magnitudes (models C, S, and V), and one case with a large coronal loop (model Y). The simulation boxes are cubes with the size of 20 Mm in “small” models (C, S and V) and 80 Mm in the “large” model Y; see Table 1 for details. The numerical grid consists of 256 × 256 × 512 elements with uniform step in each direction (i.e. the resolution in z-direction is twice as high as the resolution in the horizontal directions).

Fig. 1 Plasma temperature T and numeric density n as a function of height during the preparatory thermal relaxation phase at two different positions above the surface. Solid lines correspond to the initial state. Dotted, dashed, and dot-dashed lines represent the temporal evolution of these quantities during the thermal relaxation phase at equally spaced time intervals after the initial state. The two top panels correspond to the n and T profiles along a vertical line right above the loop footpoints, where the magnetic field is mostly vertical. The two bottom panels show the same quantities along a vertical line that crosses the centre of the loop, where the magnetic field is mostly horizontal.

All the models undergo an initial preparation phase (consisting of a thermal relaxation stage and a footpoint twisting stage). The initial magnetic field geometry is that of a potential field arcade in a gravitationally stratified atmosphere. As in Gordovskyy et al. (2014), the initial atmosphere is constructed of two isothermal regions: a lower region with thickness of about 2 Mm representing the chromosphere with temperature ≲10⁴ K and density exponentially decreasing from 5 × 10¹³ to ~10⁹ cm^-3, and an overlaying corona with temperature ~0.8 MK and nearly constant density ~10⁹ cm^-3. This idealised initial configuration is not in thermal equilibrium, however, because the thermal conduction is aligned with the magnetic field. We first let the system relax until it reaches quasi thermal equilibrium. Figure 1 shows how the vertical profiles of density and temperature evolve during the relaxation process. The temperature transition around z = 2−4 Mm becomes gradually smoother. The temperature distribution in the vertical direction tends approximately to $\mathit{T}{}^{\mathrm{\right(}}\mathit{z}^{\mathrm{\left)}}\mathrm{~}{{}^{\mathrm{(}}\mathit{z}\mathrm{+}\mathrm{const}{\mathit{.}}^{\mathrm{)}}}^{\mathrm{2}\mathit{/}\mathrm{7}}$ where the magnetic field lines are vertical. The equilibrium temperature distribution will, however, be considerably different where the field is nearly horizontal (as in the middle of the domain, across the centre of the loop arcades). Heat transfer during the relaxation phase results in some evaporation: plasma from nearly vertical field in the transition region (z = 2−3 Mm) flows upwards, noticeably increasing the loop plasma density up to z = 10 Mm. After about 2000t₀, temperature and density distributions become nearly stationary (the resulting atmosphere is thermally relaxed).

We applied footpoint twisting to the initially untwisted potential field to build the unstable twisted loop configurations. The twisting phase does not represent an actual rotation of photospheric footpoints (an effect which is often observed, but is beyond scope of this study). We can use here rotation velocities considerably higher than those observed in the photosphere to speed up the simulations, as long as these remain substantially lower than the local Alfvén speed. Under these conditions, the twisting phase should not have noticeable dynamical effects on the kink instability and reconnection (see, for more details, Gordovskyy et al. 2014, Bareford et al., 2016). The footpoint centres are located at [0.0, ± 6.5 Mm,0.0] in models C, S, and V, and at [0.0, ± 26 Mm,0.0] in model Y. Footpoints are rotated with the following profile $\mathit{\omega }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{\omega }}_{\mathrm{0}}\frac{\mathrm{1}\mathrm{-}\tanh \mathrm{[}\mathrm{-}\mathrm{(}\mathit{r}\mathrm{-}{\mathit{R}}_{\mathrm{fp}}\mathrm{)}\mathit{/}{\mathit{L}}_{\mathrm{wall}}\mathrm{]}}{\mathrm{2}}\mathit{,}$(8)where r is the distance from a footpoint centre, ω₀ ≈ 0.02 s^-1 is the frequency of rotation, R_fp is the footpoint radius, which is 0.7 Mm in models C, S, and V, and 2.8 Mm in model Y, and L_wall = 0.13 Mm is a scaling factor. Taking model S as a reference case, the footpoint rotation period is about 300 s, with a corresponding velocity at the edge of the footpoints of 1.4 × 10⁴ m/s. The coronal Alfvén speed is close to 7 × 10⁵ m/s, and the Alfvén travel time along the whole loop is ~35 s. Hence, footpoint twisting motions do not generate transients that propagate into the corona, and the coronal field has ample time to adjust and remain in near-equilibrium as it is being twisted. The other models (V and Y) have stronger magnetic fields and higher corresponding Alfvén speeds, for which quasi-steadiness is even better achieved. Chromospheric Alfvén speed is, of course, much lower, and approaches the footpoint rotation speed. However, the chromosphere has a plasma β ≫ 1, with the sound speed remaining considerably higher than the twisting speed throughout the whole chromospheric layer (up to one order of magnitude higher for model S). These reasons altogether ensure that the system evolves through a sequence of equilibria (even though small amplitude dynamical features, such as waves, can be seen during the driving period). But most importantly, the forcing (twisting) timescale is always higher than the kink instability growth-time.

Figure 5 shows the streamlines of the rotating motions applied at the surface together with the vertical magnetic field distribution B_z. These circular vortical motions are applied at the zones of maximal B_z (the footpoints), and thus are spatially coincident with the contours of the vertical field, which are approximately circular there (cf. discussions by Török & Kliem 2003; Aulanier et al. 2005).

This footpoint rotation effectively results in continued helicity injection with the total twist increasing by ~0.015 π/ s, which means that after about 500 s the total twist is of about 8π. At this stage, the twisted loop arcade is slightly thicker and taller than initially (due to the increase in magnetic pressure within the magnetic loop), but its large-scale structure stills maintains its quasi-circular symmetry. The system is close to the onset of the kink instability (but still in mechanical equilibrium), and we define this state as the initial condition (t = 0) hereafter. The configuration of the magnetic flux rope, the plasma density and temperature at this stage are represented in Fig. 2, and the corresponding pre-instability parameters are listed in Table 1.

Fig. 2 Temperature distribution shortly before the onset of the kink instability. The upper panel a) shows the temperature distribution in the mid-plane x = 0, while the lower panel b) shows the difference between the temperature distribution shortly before the instability and the initial temperature distribution. Selected magnetic field lines are shown in green in panel a). Light blue lines in panel a) denote cross-sections at which the temperature and density distributions are shown in Figs. 6 and 7. The black dashed lines in panel b) denote ΔT = 0 levels.

2.2. Synthetic soft X-ray emission from a model loop

We calculate the thermal X-ray continuum emission from the model number density n and temperature T at each point of the domain, as in Pinto et al. (2015). The plasma emissivity at a given photon energy hν is defined as $\mathit{ϵ}\mathrm{(}\mathit{h\nu ,}{\mathit{T}}^{\mathrm{)}}\mathrm{=}{\mathit{ϵ}}_{\mathrm{0}}{\mathit{n}}^{\mathrm{2}}{\mathit{T}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{g}}_{\mathrm{ff}}\mathrm{(}\mathit{h\nu ,}{\mathit{T}}^{\mathrm{)}}\exp \left(\mathrm{-}\frac{\mathit{h\nu }}{{\mathit{k}}_{\mathit{b}}\mathit{T}}\right)\mathit{,}$(9)where g_ff⁽hν,T⁾ is the Gaunt factor for free-free bremsstrahlung emission. The coefficient ϵ₀ is 6.8 × 10^-38 if the emissivity is to be expressed in units of erg cm^-3 s^-1 Hz^-1 (Tucker 1975). The Gaunt factor is approximated by the piece-wise approximation ${\mathit{g}}_{\mathrm{ff}}\mathrm{(}\mathit{h\nu ,}{\mathit{T}}^{\mathrm{)}}\mathrm{=}\{\begin{array}{c}\\ \mathrm{1}\mathit{,}& \mathit{h\nu }\lesssim {\mathit{k}}_{\mathit{b}}\mathit{T}\\ {\mathrm{(}{\frac{{\mathit{k}}_{\mathit{b}}\mathit{T}}{\mathit{h\nu }}}^{\mathrm{)}}}^{\mathrm{0.4}}\mathit{,}& \mathit{h\nu }\mathit{>}{\mathit{k}}_{\mathit{b}}\mathit{T}\mathit{.}\end{array}$(10)The photon flux density emitted at the photon energy hν is defined as $\mathit{I}\mathrm{(}\mathit{h\nu ,}{\mathit{T}}^{\mathrm{)}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}\frac{\mathrm{EM}}{\mathit{h\nu }\sqrt{{\mathit{k}}_{\mathit{b}}\mathit{T}}}{\mathit{g}}_{\mathrm{ff}}\mathrm{(}\mathit{h\nu ,}{\mathit{T}}^{\mathrm{)}}\exp \left(\mathrm{-}\frac{\mathit{h\nu }}{{\mathit{k}}_{\mathit{b}}\mathit{T}}\right)\mathit{,}$(11)where EM = n²V is the emission measure of a finite volume V of plasma of density n and temperature T. The coefficient I₀ is 1.07 × 10^-42 for a photon flux measured at a distance of 1 AU in units of photons cm^-2 s^-1 keV^-1. The total photon flux over a given spectral band is computed by integrating Eq. (11) over the corresponding range of values of hν. We compute the photon flux at different photon energies for each individual grid cell (i.e. volume element), each one having well-defined values for the emission measure, density, and temperature. As the corona is optically thin to X-ray radiation, the total flux emitted is obtained by adding the individual contributions over the whole loop (or over a region of interest).

We estimate the distributions of the differential emission measure $\mathrm{DEM}{}^{\mathrm{\right(}}\mathit{T}^{\mathrm{\left)}}$ in our simulations by computing the emission measure of the plasma regions whose temperature lies within successive temperature intervals at a given time and normalising it by the temperature interval width ΔT. This distribution then quantifies the amount of plasma at different temperatures, independently of the temperature interval considered. That is, the differential emission measure are defined here as a function of temperature as $\mathrm{DEM}\mathrm{(}{\mathit{T}}^{\mathrm{)}}\mathrm{=}\sum _{\mathit{k}}{\mathit{n}}_{\mathit{k}}^{\mathrm{2}}\mathrm{·}\frac{\mathit{\delta }{\mathit{V}}_{\mathit{k}}}{\mathrm{\Delta }{\mathit{T}}_{\mathit{k}}}\mathit{,}$(12)where the index k runs through all the plasma elements (grid cells in the simulations) that lie within the temperature interval [T,T + δT]. The values of n_k and δV_k are the number density and the volume of each element. In other words, we first compute a temperature histogram with a given temperature bin size ΔT. Then, we verify which grid cells have a temperature T within each of the bins, and we sum over all the corresponding individual EM. Variations in density for plasma at a given temperature are implicitly accounted for. Total emission measures EM can be obtained by integrating Eq. (12) relative to temperature, in the relevant interval of temperatures.

2.3. Test-particle trajectories and synthetic hard X-ray bremsstrahlung

Electron distribution functions are calculated using the GCA test-particle code based on first-order guiding-centre approximation (Gordovskyy et al. 2010; Gordovskyy & Browning 2011): $\begin{array}{ccc}& & \frac{\mathrm{d}{r}}{\mathrm{d}\mathit{t}}\mathrm{=}{u}\mathrm{+}\frac{\mathrm{\left(}\mathit{\gamma }{\mathit{v}}_{\mathrm{\left|}\mathrm{\right|}}\mathrm{\right)}}{\mathit{\gamma }}{b}\\ & & {u}\mathrm{=}{{u}}_{\mathit{E}}\mathrm{+}\frac{\mathit{m}}{\mathit{q}}\frac{\mathrm{(}\mathit{\gamma }{\mathit{v}}_{\mathrm{\left|}\mathrm{\right|}}{\mathrm{)}}^{\mathrm{2}}}{\mathit{\gamma }{\mathit{\kappa }}^{\mathrm{2}}\mathit{B}}\mathrm{[}{b}\mathrm{\times }\mathrm{(}{b}\mathrm{·}\mathrm{\nabla }\mathrm{\left)}{b}\mathrm{\right]}\\ & & \mathrm{+}\frac{\mathit{m}}{\mathit{q}}\frac{\mathit{\mu }}{\mathit{\gamma }{\mathit{\kappa }}^{\mathrm{2}}\mathit{B}}\mathrm{[}{b}\mathrm{\times }\mathrm{(}\mathrm{\nabla }\mathrm{\left(}\mathit{\kappa B}\mathrm{\right)}\mathrm{\left)}\mathrm{\right]}\\ & & \frac{\mathrm{d}\mathrm{\left(}\mathit{\gamma }{\mathit{v}}_{\mathrm{\left|}\mathrm{\right|}}\mathrm{\right)}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathit{q}}{\mathit{m}}{E}\mathrm{·}{b}\mathrm{-}\frac{\mathit{\mu }}{\mathit{\gamma }}\mathrm{(}{b}\mathrm{·}\mathrm{\nabla }\mathrm{(}\mathit{\kappa B}\mathrm{\left)}\mathrm{\right)}\\ & & \mathrm{+}\mathrm{\left(}\mathit{\gamma }{\mathit{v}}_{\mathrm{\left|}\mathrm{\right|}}\mathrm{\right)}{{u}}_{\mathit{E}}\mathrm{·}\mathrm{\left(}\mathrm{\right(}{b}\mathrm{·}\mathrm{\nabla }\mathrm{\left)}{b}\mathrm{\right)}\mathrm{+}{\left[\frac{\mathrm{d}{\mathit{v}}_{\mathrm{\left|}\mathrm{\right|}}}{\mathrm{d}\mathit{t}}\right]}_{\mathrm{coll}}\\ & & \mathit{\gamma }\mathrm{=}\sqrt{\frac{{\mathit{c}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{\gamma }{\mathit{v}}_{\mathrm{\left|}\mathrm{\right|}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{\mu B}}{{\mathit{c}}^{\mathrm{2}}\mathrm{-}{\mathit{u}}^{\mathrm{2}}}}\\ & & \frac{\mathrm{d}\mathit{\mu }}{\mathrm{d}\mathit{t}}\mathrm{=}{\left[\frac{\mathrm{d}\mathit{\mu }}{\mathrm{d}\mathit{t}}\right]}_{\mathrm{coll}}\mathrm{·}\end{array}$Here r, u, and v_|| are the particle gyro-centre position, transversal drift velocity (⊥ B) and longitudinal velocity (| | B), respectively; ${b}\mathrm{=}\frac{{B}}{\mathit{B}}$ is the magnetic field direction vector, ${{u}}_{\mathit{E}}\mathrm{=}\frac{{E}\mathrm{\times }{b}}{\mathit{B}}$ is the E × B drift component; μ is the magnetic moment per mass unit $\mathit{\mu }\mathrm{=}\frac{{\mathit{\gamma }}^{\mathrm{2}}{\mathit{u}}_{\mathrm{g}}^{\mathrm{2}}}{\mathrm{2}\mathit{B}}$, where u_g is the particle gyration velocity. The relativistic coefficients are defined as γ = (1−v²/c²)^{− 1/2} and $\mathit{\kappa }\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}{\mathit{u}}_{\mathit{E}}^{\mathrm{2}}\mathit{/}{\mathit{c}}^{\mathrm{2}}}$, where the absolute particle velocity is $\mathit{v}\mathrm{=}\sqrt{{\mathit{v}}_{\mathrm{\left|}\mathrm{\right|}}^{\mathrm{2}}\mathrm{+}{\mathit{u}}^{\mathrm{2}}\mathrm{+}{\mathit{u}}_{\mathrm{g}}^{\mathrm{2}}}$. Coulomb collisions (introduced through ${\left[\frac{\mathrm{d}{\mathit{v}}_{\mathrm{\left|}\mathrm{\right|}}}{\mathrm{d}\mathit{t}}\right]}_{\mathrm{coll}}$ and ${{}^{\mathrm{\right[}}\frac{\mathrm{d}\mathit{\mu }}{\mathrm{d}\mathit{t}}^{\mathrm{\left]}}}_{\mathrm{coll}}$) are taken into account through averaged deceleration and pitch-angle deflection of particles with mass m and charge ± e in fully ionised hydrogen plasma (see e.g. Emslie 1978). Thus,

$\frac{\mathrm{d}\mathit{v}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\mathit{a}\frac{\mathrm{2}\mathit{n}}{{\mathit{m}}^{\mathrm{2}}{\mathit{v}}^{\mathrm{2}}}\mathit{,}$with

$\mathit{a}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }{\mathit{e}}^{\mathrm{4}}}{\mathrm{(}\mathrm{4}\mathit{\pi }{\mathit{\epsilon }}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}\mathrm{\Lambda }\left(\frac{{\mathit{m}}_{\mathrm{e}}\mathrm{+}{\mathit{m}}_{\mathrm{p}}}{{\mathit{m}}_{\mathrm{e}}{\mathit{m}}_{\mathrm{p}}}\right)\mathit{,}$where Λ is the Coulomb logarithm, and quasi-random changes of particle pitch-angles with probabilities, and, hence, $\mathit{v}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}\mathrm{\Delta }\mathit{v}\mathrm{=}\mathrm{3}\sqrt{{\mathit{v}}^{\mathrm{3}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{-}\frac{\mathrm{2}\mathit{an}}{{\mathit{m}}^{\mathrm{2}}}\mathrm{\Delta }\mathit{t}}\mathit{.}$(20)Corresponding quasi-random angular deflections are calculated as $\mathrm{\Delta }\mathit{\theta }\mathrm{=}\mathrm{1.3}\frac{\mathit{an}}{{\mathit{m}}^{\mathrm{2}}{\mathit{v}}^{\mathrm{3}}}\sqrt{\mathit{\mu }}\mathrm{𝒩}\sqrt{\mathrm{\Delta }\mathit{t}}\mathit{,}$(21)where is a random number from normal distribution with the width of 1.0. Tests show that, statistically, this formula provides a good approximation for pitch-angle diffusion as per $\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{=}\frac{\mathit{an}}{\mathit{m}{\mathit{v}}^{\mathrm{3}}}\frac{\mathit{\partial }}{\mathit{\partial \mu }}\left(\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}\frac{\mathit{\partial }}{\mathit{\partial \mu }}\right)$, where μ = cosθ = v_||/v is the pitch-angle cosine (Landau 1937). Therefore, collisional terms can be written using Eqs. (20), (21): $\begin{array}{ccc}{\left[\frac{\mathrm{d}\mathit{\mu }}{\mathrm{d}\mathit{t}}\right]}_{\mathrm{coll}}& \mathrm{=}& \sqrt{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}\frac{\mathrm{\Delta }\mathit{\theta }}{\mathrm{\Delta }\mathit{t}}\\ {\left[\frac{\mathrm{d}{\mathit{v}}_{\mathrm{\left|}\mathrm{\right|}}}{\mathrm{d}\mathit{t}}\right]}_{\mathrm{coll}}& \mathrm{=}& \mathit{\mu }\frac{\mathrm{\Delta }\mathit{v}}{\mathrm{\Delta }\mathit{t}}\mathrm{+}\mathit{v}{\left[\frac{\mathrm{d}\mathit{\mu }}{\mathrm{d}\mathit{t}}\right]}_{\mathrm{coll}}\mathrm{·}\end{array}$Initially, the box is randomly filled with 2²⁰ test-electrons that are uniformly distributed in space and relative to pitch-angle cosine μ, while relative to the absolute velocity v particles have a Maxwellian distribution corresponding to T = 0.83 MK (the temperature of the coronal background; see Fig. 1). Thermal bath conditions are applied on all six boundaries of the domain: particles are free to leave the domain, and for every particle leaving the domain, one new particle is generated randomly from a Maxwellian distribution (at the thermal bath temperature) and is injected at the same point of the boundary.

Hard X-ray bremsstrahlung intensities are calculated using the Kramers formula (see e.g. Kontar et al. 2002),

$\mathit{I}\mathrm{\left(}{\mathit{ϵ}}_{\mathrm{ph}}\mathrm{\right)}\mathrm{=}{\mathit{Q}}_{\mathrm{0}}\underset{{\mathit{ϵ}}_{\mathrm{ph}}}{\overset{\mathrm{\infty }}{\mathrm{\int }}}{\mathit{f}}_{\mathrm{beam}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathit{n}\frac{\mathrm{1}}{{\mathit{ϵ}}_{\mathrm{ph}}\mathit{ϵ}}\mathrm{d}\mathit{ϵ,}$where I(ϵ_ph) is the rate of emission of photons with energy ϵ_ph, ϵ is an electron energy, f_beam is the differential beam electron density (dn/ dE), n is the space and time-varying electronic density of the ambient plasma, and Q₀ is the constant from the Kramers bremsstrahlung cross-section formula.

The test particle approach is widely used to study particle acceleration in flares, but is not a fully self-consistent approach since it neglects the electromagnetic fields generated by the test particles themselves. However, in the configurations studied here, the fraction of the released magnetic energy carried by the accelerated particles is relatively small – typically around 5% (Gordovskyy et al. 2014). Furthermore, as a result of the fragmented nature of the current sheet, energetic electrons are fairly uniformly distributed within the loop volume, so we expect that their effects on the fields are weak. Hence, the test particle approach is quite well justified.

Fig. 3 Magnetic field topology and current distribution at four representative instants of the simulations (from top to bottom: t = 0, 30, 120, 196 s after the onset of the kink instability; model S), illustrating the onset and growth of the kink instability, the reconnection episode that follows, and the relaxation of the coronal loop towards a non-twisted state. The lines in the left column are magnetic field lines, with the blue lines having seed points very close to the centres of the flux rope footpoints. The axes are labelled in units of Mm. The right column shows iso-surfaces of the current density j = jcrit (see Eq. (7) and the following discussion).

Fig. 4 Evolution of the magnetic, internal, and kinetic energy in model S after the onset of the kink instability.

Fig. 5 Footpoint twisting profile and vertical magnetic field amplitude at the lower boundary (for model S). The grey-scale shows the amplitude of the vertical magnetic field Bz (white and black representing positive and negative polarity) and the blue lines are streamlines of the surface motions (see Eq. (8)). The axes are labelled in units of Mm.

Fig. 6 Temperature distributions during and after the kink instability. Three loop cross-sections are shown (corresponding to the blue lines in Fig. 2): left footpoint (z = 1.9 Mm) (left column), left loop leg (z = 4.25 Mm) (middle column), and vertical plain crossing the loop top (right column). Times shown are (after onset of the kink instability): panel a) – 0 s, b) – 81 s, c) – 196 s, and d) – 377 s.

Fig. 7 Logarithm of density distributions. Locations and times are the same as in Fig. 6.

3. Magnetic fields, currents, and hot plasma

The kink instability in the present MHD simulations develops when the total twist reaches around 7−9π. It leads to the formation of a helical current sheet around the twisted flux-rope (see Fig. 3), in a way similar to that described by Brown et al. (2003), Browning et al. (2008), Hood et al. (2009), Gordovskyy & Browning (2012), Gordovskyy et al. (2014), and Pinto et al. (2015). This in turn switches the anomalous resistivity on (where and when the current density becomes stronger than j_crit; see Eq. (7)) and leads to the start of the reconnection in the corona. There also are very strong currents near the loop footpoints, but the critical current threshold j_crit in the lower atmosphere is much higher, and hence the resistivity there remains equal to the much lower background resistivity. Hence, the reconnection and current dissipation occur predominantly in the coronal part of a loop. The current density distribution evolves gradually into a more fragmented pattern thereafter, until it fades away. Magnetic reconnection occurs both within the twisted loop and between twisted flux rope field lines and the ambient field.

The reconnection episode is followed by a magnetic relaxation phase lasting about 100−200 s in the smaller models and about 300−400 s in the large-scale models. During this time, the magnetic energy is reduced by ~50−70% of the free magnetic energy in the system (Fig. 4). The amount of free magnetic energy is defined here as the difference between the magnetic energy just before kink instability and the initial magnetic energy of the potential (untwisted) magnetic configuration. Most of this energy (around 80−90%) is converted directly into thermal energy through ohmic dissipation, while a small fraction is released in form of plasma kinetic energy. The latter is then also converted into thermal energy through viscous effects, mostly around shocks. A more detailed description of the velocity field in reconnecting twisted loops is given by Gordovskyy et al. (2015).

Ohmic dissipation is much enhanced in the strong helical current, leading to plasma heating. Figures 6 and 7 show the spatial distribution of the plasma temperature and density at different stages of the evolution of the system. Initially, the plasma temperature and density are everywhere equal to those of the initial gravitationally-stratified background (cf. Fig. 1). During the thermal relaxation stage, thermal conduction increases the temperature and decreases the density in the transition region while reducing the temperature and increasing the density in the lower corona. Footpoint twisting leads to an increase in magnetic pressure inside the twisted loop, causing it to expand. This expansion in turn results in a slight drop in plasma density and temperature inside the loops (by adiabatic expansion). Just before the kink instability occurs, and near the loop-top, the temperature is about 40−50% lower than the temperature in the initial gravitationally stratified atmosphere (see Figs. 2 and 6a), while the density is 50−60% lower (Fig. 7a). The picture is more complicated closer to the foot-points: plasma around the expanding twisted flux tube is being pushed vertically and not horizontally by the strong vertical density gradient. This results in accumulation of mass in a cylindrical shell around the loop legs. The density in these shells is about 2.5 times higher than outside, while the temperature is about 1.5−2 times higher because of adiabatic compression. After the reconnection starts, the high temperature distributions near the loop top resemble the current density distributions (e.g. Hood et al. 2009), starting with a helical shell, which then gradually becomes more fragmented and fades away. In contrast, closer to foot-points the hot plasma can be seen on the central axis of the loop. The density near the loop top evolves from a regular structure towards a fragmented one. In addition, a substantial increase in plasma density after onset of reconnection can be observed in the loop legs, which is due to plasma evaporation from the lower atmospheric layers.

Figure 8 shows the $\mathrm{DEM}{}^{\mathrm{\right(}}\mathit{T}^{\mathrm{\left)}}$ of the coronal plasma at different instants of the simulation for three different models (S, V, and Y). The DEM distribution extends into higher temperatures during the linear phase of the kink instability in all the models, forming a distinctive upper tail. The total emission measure of this hotter plasma component continues to increase until the peak phase is reached, and a clear hot peak is sometimes visible (especially in model V). The bulk of this hotter plasma reaches temperatures of about twice the temperature of the background plasma in these simulations. The DEM distribution flattens out afterwards, during the magnetic relaxation phase, but the maximum temperature remains higher than initially for a long period of time. Overall, the temporal evolution of the DEM analysed here is very similar to those discussed by Pinto et al. (2015) for a simpler model of a kink-unstable flaring loop. They related the formation and the dissipation of hot flare plasma components in the DEM to observations by Reale et al. (2009), Battaglia & Kontar (2012) and Sylwester et al. (2014), for example. The only remarkable differences are that in our case the flare plasma reaches lower maximum temperatures and that the hot plasma component that forms during the course of the flare is less pronounced (because the ratio of background to flare plasma is higher here as well).

4. SXR and HXR emission

The combination of the variations in plasma density and temperature gives rise to the soft X-ray thermal emission patterns displayed in Figs. 9. The figure shows volume renderings of the thermal continuum emissivity at 2 keV at several instants (t = 34, 80, 136, and 227 s, respectively, from top to bottom) together with renderings of the plasma temperature. The left and right columns of each subplot show the same loop from different points of view (side and perspective, and top views). The orange and blue lines represent magnetic field lines rooted at the left and right foootpoints. In the left subplot, the yellow to orange volumes represent the plasma emissivity ϵ at 2 keV (see Eq. (9)), ranging from ϵ = 5 × 10¹² (light yellow) to ϵ = 5 × 10¹⁴ erg s^-1 cm^-3 Hz^-1 (dark orange). In the right subplot, the light yellow to dark orange volumes represent temperatures ranging from 1 MK to 2 MK. The continuum thermal plasma emissivity does not exactly follow the evolution of the plasma temperature. The difference is due to the variations in plasma density that result from the turbulent dynamics of the kinking loop – we note the dependence on density in Eq. (9) –, as was also the case in Pinto et al. (2015).

The first traces of thermal emission appear at mid-height and around the loop legs. These occur in the form of two separate sheet-like emission structures placed around the boundaries of the loop, which are pushed against the background medium as the flux ropes writhes. Thermal emission then fills up the whole flux rope volume for a short period of time (at about the maximum of emission). The total thermal emission fades out afterwards as the loop plasma cools down globally during the relaxation phase. At this point, the thermal emission concentrates into a few separate field-aligned threads, most notably at mid-heights and crossing the top part of the loop (see, especially, the top views in Fig. 9). As the figure shows, thermal emission is only clearly discernible well past the onset of the kink instability, when the flux rope has already lost an important fraction of its initial twist. Therefore, the amount of twist that can be measured by tracing out the emission threads is much lower than the value of the initial magnetic twist. The two quantities differ by a factor 4 in our simulations (the emission threads show a twist of about 2π, while the initial magnetic twist is ~8π). This situation is analogous to that described by Pinto et al. (2015) for a simplified model of a kink-unstable coronal loop (without large-scale curvature). We note that the thermal emission shows a structure with some writhe, which is in the sense of the original field line twist. Figure 10 represents the thermal photon flux at the same photon energy (2 keV) and at the same times as Fig. 9, but integrated along the line of sight (l.o.s.) for three different viewpoints (side, front, and top of the loop), with the darker tones representing stronger photon flux. The contour lines show the locations of the peaks of non-thermal hard X-ray bremsstrahlung (HXR) emission at 5 and 12 keV (blue and green contours, respectively, at 90% and 75% of the peak emission). The hard X-ray emission spatial profiles are smoothed with a Gaussian filter with ~1 Mm width, roughly corresponding to the resolution of modern instruments (we did not try to reproduce a particular instrument performance, however). It is interesting to note that depending on the l.o.s., the thermal emission may appear as coming from a unique coronal loop or from two independent sources. This occurs because the emission structures that are aligned with the l.o.s. are much enhanced with respect to the emission structures oriented orthogonally because the former occupies a shorter depth along the l.o.s. than the latter. The top views (rightmost columns) in Figs. 9 and 10 show this effect very clearly. The emission threads that cross the loop above its apex are well visible in the volume renderings in Fig. 9, but nearly disappear when integrated along a vertical l.o.s. (Fig. 10). This effect further hinders the possibility of measuring the actual levels of magnetic field line twist directly from the morphology of the thermal emission threads.

Fig. 8 Temporal evolution of the DEM for cases S, V and Y (left, middle and right columns, respectively). Time runs from the top to the bottom rows (as indicated in the inset legend on each individual panel, where time is indicated in units of t0). For each one of the cases represented, the first row corresponds to the linear phase of the kink-instability, the second row represents the period of maximum thermal emission, and the third row the beginning of the relaxation phase. The upper tail which forms in the DEM distribution corresponds to the plasma heated after the development of the kink-instability in the magnetic flux rope. This corresponds to a hot plasma component, which is particularly distinguishable in case V. The maximum temperature reached depends mostly on flux rope’s magnetic field amplitude, approximately as ~ B2 (cf. Pinto et al. 2015).

Fig. 9 Volume renderings of the thermal continuum emissivity at 2 keV (left) and of the plasma temperature (right) at different instants for model S. The left and right columns of each plot show a perspective view and a top view of the same loop. The colour scales represent emissivities (ϵ; see Eq. (9)) between 5 × 1012 and 5 × 1014 erg s-1 cm-3 Hz-1 and plasma temperatures between 1 and 2 MK. The orange and blue lines represent magnetic field lines rooted at the left and right footpoints. The instants represented are from top to bottom, t = 34, 80, 136 and 227 s. The grey line delimits the bottom face of the numerical domain (a square of dimensions 20 × 20 Mm).

Fig. 10 Thermal continuum photon flux emitted at 2 keV (colour-scale) and hard X-ray emission at 5 and 12 keV (blue and green contours, respectively, at 90% and 75% of the peak emission at those photon energies) integrated along the l.o.s. in three different directions (from left to right: side, edge, and top views) for our reference case (model S). The thermal and non-thermal photon fluxes are normalised to arbitrary units to facilitate comparison. The axes are in units of Mm. The instants represented are from top to bottom, t = 34, 80, 136, and 227 s (as in Fig. 9).

Fig. 11 Spatial distribution of energetic electrons in the model S seen from the loop side (cf. the first column of Fig. 10) at different instants. The left a), middle b), and right c) columns correspond to 35 s, 81 s, 104 s after onset of the kink instability. The axes indicate distances in Mm. Red and green colours indicate high and low electron densities, respectively. Upper and lower panels correspond to particles with energies > 5 keV and > 12 keV, respectively.

Fig. 12 Top row: light curves of the soft X-ray emission (black and red lines) and hard X-ray emission (blue and green lines) in models C, S, and V at different photon energy bands (see the inset legend). The bottom row shows the temporal derivative of the soft X-ray light curves together with the hard X-ray light curves during the first 150 s for the same models. The curves are all normalised to their peak value to facilitate representation. The peaks in the higher energy SXR light curves (red lines) tend to precede the lower energy ones (black lines); this is especially visible in models C and S. The decay is also faster at higher SXR photon energy bands. The HXR light curves (blue and green lines) tend to peak at the moment when the SXR growth rate is at its highest (this is especially visible in model S), such that the temporal derivatives of the SXR curves correlate well with the corresponding HXR curves. A few post-flare oscillations are visible in all the SXR light curves represented.

Generally, the non-thermal electron behaviour we observe here is similar to that described in Gordovskyy et al. (2014). However, electrons appear to lose their energies faster. This most likely happens because of the higher density closer to the footpoints. Thus, thermal conduction in the initial atmosphere in the present model leads to an increase of the density by factor of 2 around z ≈ 3−5 Mm over footpoints and by factor of ~3 around z ≈ 2−5 Mm in horizontal field (see Fig. 1). In addition, dense plasma is also accumulated around loop legs during the twisting (Fig. 7). Similarly, hard X-ray bremsstrahlung emission evolves slightly faster than in previous simulations. Although there are very many energetic electrons throughout the coronal loop, the bulk of the emission comes from very flat footpoint sources with a height of 1−2 Mm (comparable to the resolution of these intensity maps). But the situation might in principle be different in coronal loop configurations with higher coronal densities (HXR emission peaks might eventually appear at higher coronal heights, as is sometimes observed in solar flares). There HXR footpoint emission occasionally shows slight asymmetries, but these disappear quite quickly. Overall, the HXR emission is essentially similar and simultaneous in both footpoints. The area of HXR footpoints initially increases with time as a result of the twisted loop reconnection with ambient field (see Gordovskyy & Browning 2012). This increase is not monotonic, and at some point, the footpoint areas begin to shrink. This is the consequence of two different factors: numerical undersampling due to the limited test-particle number, and modulation of the HXR intensities. Figure 11 show the spatial distribution of electrons with energies above 5 and 12 keV at different moments of the temporal evolution of model S (the same as in Fig. 10). Electron distribution functions are calculated by reducing test-particles on a 40 × 40 × 40 × 40[x,y,z,ϵ] grid with uniform spatial step and exponential energy step to reduce undersampling at higher energies. This is clearly very similar to the distribution of the parallel electric field in the domain (see Fig. 3), and hence energetic electrons are distributed within most of the flaring loop volume, even though the HXR emission in Fig. 10 is much more concentrated in the lower layers. We note that the HXR emission spatial patterns are not a direct reflection of the electron acceleration mechanisms, but rather of the interaction between the accelerated electrons with the highly stratified coronal plasma. A substantial amount of energetic particles accumulates, nonetheless, near the loop footpoints. This is because the magnetic field convergence reduces the parallel component of the particle velocities. The magnetic trapping is even more effective because of collisional pitch-angle scattering that moves particles out of the loss-cone. Since all the magnetic field lines of the twisted loop are tied to the lower boundary (i.e. the photosphere), virtually all energetic particles leaving the domain do so through the lower boundary.

Figure 12 shows the light curves of the total thermal continuum emission integrated in three photon energy bands between 1 and 12 keV (see inset legend) and of the hard X-ray emission at 5 and 12 keV for models C, S, and V. The light curves are all normalised to their peak value for clarity (absolute photon fluxes decrease with increasing photon energy). The emission in lower energy bands is sustained for longer times (compared to non-thermal emission), as is commonly observed in flares. Photon fluxes at lower energy bands (e.g. 1−3 keV) systematically have a slow rise and decay, while flux at higher energy bands (e.g. 6−12 keV) always displays a more rapid evolution (as in Pinto et al. 2015). The peak photon fluxes do not occur simultaneously at all photon energy bands: the peaks in higher energy light curves tend to precede those in lower energy light curves. This is particularly noticeable in models S and V. A few post-flare oscillations are visible in all light curves represented (hence, in an energy range between 1 and 25 keV). As expected, non-thermal emission evolves much faster than thermal X-ray emission. In model S, for example, the HXR flux peaks at about t = 60 s, already after the saturation of the kink instability and the triggering of the fast reconnection, when the initial large-scale helical current sheet has already started fragmenting. This is visible by comparison with Fig. 3; the HXR peak occurs just in between the second and the third panels (t = 30 and 120 s, respectively). The current density distribution is indeed very much fragmented during this whole time interval, even though its large-scale spatial organisation is a remnant of the initial helical current sheet. Both 5 and 12 keV bremsstrahlung emission decay thereafter and nearly disappear after ~100−200 s, when thermal emission is just reaching its peak values.

5. Discussion

We study here the consequences of the triggering of the kink instability in twisted coronal loops in terms of X-ray emission (thermal continuum SXR and bremsstrahlung HXR) by means of numerical MHD simulations combined with test-particle methods. The MHD model consists of a semi-circular kink-unstable coronal loop embedded in a stratified atmosphere composed of coronal and chromospheric layers. The resulting plasma density and temperature are used to calculate emission measures and thermal continuum emission (soft X-rays), and the electric and magnetic fields are used to compute test-particle trajectories. Non-thermal (hard X-ray) emission is estimated based on the interaction between these test-particles and the stratified background. Heating caused by these collisions and currents generated by the accelerated particles are not fed back into the MHD simulation. That would require a different type of approach based on hybrid methods, which is clearly beyond the scope of this manuscript.

The vertical grid resolution – ~0.16 Mm and ~0.04 Mm for the “large” and “small” models, respectively – under-resolves the steep transition region (TR) between the chromosphere and the corona, even though they are at the forefront of current simulations of the solar corona. For comparison, our vertical resolution is similar to or higher than the landmark three-dimensional simulations by Gudiksen & Nordlund (2005) and Bingert & Peter (2011) at 0.15 Mm and 0.23 Mm, respectively, and approach the vertical resolution requirement at the TR of 20 km suggested by Klimchuk (2006). Under-resolving the TR can have an impact on the properties of the simulated corona (density, temperature, and/or heat fluxes), hence possibly affecting both the synthetic SXR and HXR emissions. Bradshaw & Cargill (2013) have shown that under-resolving the transition region may significantly underestimate the peak density of coronal loops undergoing impulsive heating episodes (and being loaded with plasma by means of chromospheric evaporation), but without much effect on the plasma temperatures they reach. In spite of these limitations, our results are consistent with the main common features of observed solar flares and furthermore allow us to clearly establish links between different types of X-ray emission.

Our models represent confined (non-eruptive) flares in single loops. The cases studied here represent smaller – and also thus more frequent – events. Larger energy releases can be obtained by considering stronger magnetic fields and larger loops (Gordovskyy et al. 2015). However, whilst the kink instability in a loop may play a role in the largest flares (e.g. by triggering a filament eruption), such events would involve larger-scale and more complex field configurations than considered here. Our results provide a prediction of what is expected to be observed in kink-unstable loops, and thus new observational campaigns are required to look for such signatures. Detailed comparison between simulations and observations is a task for future work. However, much work already published provides evidence of events that appear to be consistent with our models. The simulated loops rise somewhat during the main reconnection phase, then contract back down, so may be observed as “failed eruptions”, which appear sometimes to be associated with kink instability (e.g. Alexander et al. 2006; Song et al. 2014). Flares can be directly associated with the apparent instability of observed twisted flux ropes (e.g. Wang et al. 2015). Liu et al. (2013) analysed a flare within a loop, suggesting that the x-ray emission is generated by fragmented current sheets in the coronal part of the loop. Furthermore, there is evidence for an association of a period of helicity build-up before flare onset (Reinard et al. 2010), consistent with our scenario.

Our results suggest that highly twisted coronal loops may be more abundant than current observations seem to indicate. Individual magnetic field lines in flaring coronal loops will only be visible as emission threads if (and only if) the enclosed plasma has already been considerably heated with respect to the background plasma during the course of the flare. At this point, magnetic field lines twisted well above the threshold for the kink instability will already have lost an important fraction of their initial twist. The effect is here even stronger than in the simulations reported by Pinto et al. (2015), in that part of the emission threads nearly disappear in the l.o.s. projections. Comparing Fig. 10, which shows the photon fluxes integrated along the l.o.s. for the three-dimensional emission structures shown in Fig. 9, makes this particularly evident. The threads that cross the loop above its apex nearly disappear when integrated along a vertical line of sight, for example. These emission patterns could even be interpreted as coming from a coronal structure with a slight large-scale torsion but without any twist at all. This result removes one of the main obstacles faced by the kink-instability scenario for the theoretical understanding of the triggering and evolution of solar flares.

The temporal ordering of the soft and hard X-ray light curves we obtained are consistent with those typically observed. The light curves in Fig. 12 show similarities with GOES and RHESSI light curves, such as those reported by Sylwester et al. (2014) for an M1.0 class flare. Lower energy thermal (SXR) light curves grow and decay more slowly than higher energy thermal light curves (Pinto et al. 2015, in a way similar to that in). Furthermore, the hard X-ray (HXR) light curves that are due to the bremsstrahlung of non-thermal electrons reach their peak value when the growth-rate of the SXR curves are maximal. This behaviour is consistent across all the models we tested (C, S, V, and Y) and approximately reproduces the so-called Neupert effect. The most commonly accepted interpretation of this effect invokes heating of chromospheric plasma by the supra-thermal electron beams that produce HXR emission, followed by the up-flow of hot plasma into the coronal loops and subsequent SXR emission. However, our simulations suggest a different cause for this effect, as the electron beams we consider are not allowed to heat up the background plasma by construction (the supra-thermal electrons are test-particles). Here, the instant at which a given HXR light curve peaks is very close to the time when the current sheet that forms around the twisted coronal loop reaches its maximum amplitude and spatial extent (i.e. just before the main magnetic reconnection episode starts). Particle acceleration is maximal at this instant, and so are the flux of particle reaching the dense lower layers of the atmosphere and the consequent HXR photon flux. Then, actual plasma heating leading to SXR emission continues to occur afterwards during the relaxation phase (even though the ohmic heating is proportional to ηj², and hence maximal at that instant). This result agrees with some of the conclusions of Veronig et al. (2005), namely that fast electron beams may not be the main source of SXR plasma supply and heating, and that other heating sources are necessary to explain the Neupert effect quantitatively.

The ensemble of results presented in this manuscript support the view that plasma heating and particle acceleration processes in solar flares may occur in large volumes of the respective flaring loops instead of in a small volumes located near their apexes. The flare scenario based on the triggering of the kink instability of coronal loops with magnetic helicity (twist) naturally provides the physical mechanisms necessary for volume-distributed heating and acceleration.

6. Summary

The aim of this study was to determine observational manifestations of solar flares occurring in kink-unstable twisted coronal loops. Specifically, we focused on thermal soft X-ray (SXR) continuum emission together with non-thermal hard X-ray (HXR) bremsstrahlung emission deduced from a set of MHD simulations. The thermal SXR emission was calculated from the distributions of density and temperature of the compressible plasma, and the HXR emission was derived from the interaction between test particles and the time-evolving MHD stratified background.

Our results showed that the geometry of the emission produced during flares does not necessarily correspond to that of the underlying magnetic field. In particular, the twist angles revealed by the emission threads (soft X-ray thermal emission; SXR) are consistently lower than the actual field line twist present at the onset of the kink instability. Hard X-ray (HXR) emission due to the collisions of energetic electrons with the stratified background is concentrated at the loop footpoints (quasi-simultaneously), even though the electrons are accelerated everywhere within the coronal volume of the loop. The maximum of HXR emission consistently precedes that of SXR emission. The HXR light curve is approximately proportional to the temporal derivative of the SXR light curve, hence reproduces the so-called Neupert effect.

Future work will cover other aspects of thermal emission in flares such as thermal emission in the extreme ultraviolet energy bands (continuum and line emission), and also properties of micro-wave gyro-synchrotron emission in flares triggered by the kink instability.

Acknowledgments

We thank the anonymous referee, whose insightful remarks helped us improving this manuscript considerably. R.F.P. was funded by the French National Space Centre (CNES) and by the FP7 project #606692 (HELCATS). M.G. and P.K.B. are funded by the Science and Technology Facilities Council (UK). Numerical simulations have been performed using the UKMHD computational facilities (part of the STFC DiRAC HPC).

References

All Tables

All Figures

Fig. 1 Plasma temperature T and numeric density n as a function of height during the preparatory thermal relaxation phase at two different positions above the surface. Solid lines correspond to the initial state. Dotted, dashed, and dot-dashed lines represent the temporal evolution of these quantities during the thermal relaxation phase at equally spaced time intervals after the initial state. The two top panels correspond to the n and T profiles along a vertical line right above the loop footpoints, where the magnetic field is mostly vertical. The two bottom panels show the same quantities along a vertical line that crosses the centre of the loop, where the magnetic field is mostly horizontal.

In the text

Fig. 2 Temperature distribution shortly before the onset of the kink instability. The upper panel a) shows the temperature distribution in the mid-plane x = 0, while the lower panel b) shows the difference between the temperature distribution shortly before the instability and the initial temperature distribution. Selected magnetic field lines are shown in green in panel a). Light blue lines in panel a) denote cross-sections at which the temperature and density distributions are shown in Figs. 6 and 7. The black dashed lines in panel b) denote ΔT = 0 levels.

In the text

Fig. 3 Magnetic field topology and current distribution at four representative instants of the simulations (from top to bottom: t = 0, 30, 120, 196 s after the onset of the kink instability; model S), illustrating the onset and growth of the kink instability, the reconnection episode that follows, and the relaxation of the coronal loop towards a non-twisted state. The lines in the left column are magnetic field lines, with the blue lines having seed points very close to the centres of the flux rope footpoints. The axes are labelled in units of Mm. The right column shows iso-surfaces of the current density j = jcrit (see Eq. (7) and the following discussion).

In the text

	Fig. 4 Evolution of the magnetic, internal, and kinetic energy in model S after the onset of the kink instability.
In the text

	Fig. 5 Footpoint twisting profile and vertical magnetic field amplitude at the lower boundary (for model S). The grey-scale shows the amplitude of the vertical magnetic field Bz (white and black representing positive and negative polarity) and the blue lines are streamlines of the surface motions (see Eq. (8)). The axes are labelled in units of Mm.
In the text

Fig. 6 Temperature distributions during and after the kink instability. Three loop cross-sections are shown (corresponding to the blue lines in Fig. 2): left footpoint (z = 1.9 Mm) (left column), left loop leg (z = 4.25 Mm) (middle column), and vertical plain crossing the loop top (right column). Times shown are (after onset of the kink instability): panel a) – 0 s, b) – 81 s, c) – 196 s, and d) – 377 s.

In the text

	Fig. 7 Logarithm of density distributions. Locations and times are the same as in Fig. 6.
In the text

Fig. 8 Temporal evolution of the DEM for cases S, V and Y (left, middle and right columns, respectively). Time runs from the top to the bottom rows (as indicated in the inset legend on each individual panel, where time is indicated in units of t0). For each one of the cases represented, the first row corresponds to the linear phase of the kink-instability, the second row represents the period of maximum thermal emission, and the third row the beginning of the relaxation phase. The upper tail which forms in the DEM distribution corresponds to the plasma heated after the development of the kink-instability in the magnetic flux rope. This corresponds to a hot plasma component, which is particularly distinguishable in case V. The maximum temperature reached depends mostly on flux rope’s magnetic field amplitude, approximately as ~ B2 (cf. Pinto et al. 2015).

In the text

Fig. 9 Volume renderings of the thermal continuum emissivity at 2 keV (left) and of the plasma temperature (right) at different instants for model S. The left and right columns of each plot show a perspective view and a top view of the same loop. The colour scales represent emissivities (ϵ; see Eq. (9)) between 5 × 1012 and 5 × 1014 erg s-1 cm-3 Hz-1 and plasma temperatures between 1 and 2 MK. The orange and blue lines represent magnetic field lines rooted at the left and right footpoints. The instants represented are from top to bottom, t = 34, 80, 136 and 227 s. The grey line delimits the bottom face of the numerical domain (a square of dimensions 20 × 20 Mm).

In the text

Fig. 10 Thermal continuum photon flux emitted at 2 keV (colour-scale) and hard X-ray emission at 5 and 12 keV (blue and green contours, respectively, at 90% and 75% of the peak emission at those photon energies) integrated along the l.o.s. in three different directions (from left to right: side, edge, and top views) for our reference case (model S). The thermal and non-thermal photon fluxes are normalised to arbitrary units to facilitate comparison. The axes are in units of Mm. The instants represented are from top to bottom, t = 34, 80, 136, and 227 s (as in Fig. 9).

In the text

Fig. 11 Spatial distribution of energetic electrons in the model S seen from the loop side (cf. the first column of Fig. 10) at different instants. The left a), middle b), and right c) columns correspond to 35 s, 81 s, 104 s after onset of the kink instability. The axes indicate distances in Mm. Red and green colours indicate high and low electron densities, respectively. Upper and lower panels correspond to particles with energies > 5 keV and > 12 keV, respectively.

In the text

Fig. 12 Top row: light curves of the soft X-ray emission (black and red lines) and hard X-ray emission (blue and green lines) in models C, S, and V at different photon energy bands (see the inset legend). The bottom row shows the temporal derivative of the soft X-ray light curves together with the hard X-ray light curves during the first 150 s for the same models. The curves are all normalised to their peak value to facilitate representation. The peaks in the higher energy SXR light curves (red lines) tend to precede the lower energy ones (black lines); this is especially visible in models C and S. The decay is also faster at higher SXR photon energy bands. The HXR light curves (blue and green lines) tend to peak at the moment when the SXR growth rate is at its highest (this is especially visible in model S), such that the temporal derivatives of the SXR curves correlate well with the corresponding HXR curves. A few post-flare oscillations are visible in all the SXR light curves represented.

In the text