3D coupled tearing-thermal evolution in solar current sheets

Combined tearing-thermal evolution plays an important role in the disruption of current sheets, and formation of cool condensations within the solar atmosphere. However, this has received limited attention to date. We numerically explore a combined tearing and thermal instability that causes the break up of an idealized current sheet in the solar atmosphere. The thermal component leads to the formation of localized, cool condensations within an otherwise 3D reconnecting magnetic topology. We construct a 3D resistive magnetohydrodynamic simulation of a force-free current sheet under solar atmospheric conditions that incorporate the non-adiabatic influence of background heating, optically thin radiative energy loss, and magnetic field aligned thermal conduction with the open source code MPI-AMRVAC. Multiple levels of adaptive mesh refinement reveal the self-consistent development of finer-scale condensation structures within the evolving system. The instability in the current sheet is triggered by magnetic field perturbations concentrated around the current sheet plane, and subsequent tearing modes develop. This in turn drives thermal runaway associated with the thermal instability of the system. We find subsequent, localized cool plasma condensations that form under the prevailing low plasma-$\beta$ conditions, and demonstrate that the density and temperature of these condensed structures are similar to more quiescent coronal condensations. Synthetic counterparts at Extreme-UltraViolet (EUV) and optical wavelengths show the formation of plasmoids (in EUV), and coronal condensations similar to prominences and coronal rain blobs in the vicinity of the reconnecting sheet. Our simulations imply that 3D reconnection in solar current sheets may well present an almost unavoidable multi-thermal aspect, that forms during their coupled tearing-thermal evolution.


Introduction
Magnetic reconnection is a fundamental process understood to play a critical role throughout the solar atmosphere.The change of magnetic field topology during reconnection leads to conversion of magnetic energy into thermal and kinetic energies (Biskamp 2000), frequently leading to fast energy release of solar flares (Giovanelli 1939(Giovanelli , 1947(Giovanelli , 1948;;Priest & Forbes 2000;Hesse & Cassak 2020) and coronal mass ejections (Gosling et al. 1995;Schmidt & Cargill 2003;Karpen et al. 2012) out into the heliosphere.
It was suggested by Furth et al. (1963) that reconnection in an incompressible plasma may be triggered due to small perturbations in a current layer, correspondingly breaking up the current sheet in the form of the tearing instability.Linear analysis by Loureiro et al. (2007) suggests that the tearing instability in a single current sheet may lead to the formation of a chain of plasmoids (secondary magnetic islands).This was later verified in 2D numerical simulations by Huang & Bhattacharjee (2013) and Huang et al. (2013).Extensions to 2D double current layer models were seen to give rise to the development and layer-layer interactions of tearing modes with smaller scale plasmoid formation (Zhang & Ma 2011;Keppens et al. 2013;Akramov & Baty 2017;Paul & Vaidya 2021, and references therein).
Each of these previous studies does not consider nonadiabatic effects of background heating, radiative energy loss, and thermal conduction, which are essential components of the solar atmosphere.It is well established that the solar corona is in an overall delicate thermal balance.If this balance due to optically thin radiative loss, background heating in combination with thermal conduction is perturbed, an increment of the thermal energy loss cooling down the plasma may lead to an enhancement of the plasma density.This in turn radiates more energy (radiative loss in an optically thin medium is proportional to the square of plasma density), and the material becomes cooler still.Hence, a catastrophic runaway process ensues leading to a rapid rise in the density and drop in temperature which is in essence the thermal instability.A detailed linear analysis of the thermal instability is presented by Parker (1953) and Field (1965), who derived the criteria governing the onset to a catastrophic radiative loss in an infinite homogeneous medium.The linear magnetohydrodynamic (MHD) analysis was extended to a 1D slab configuration (van der Linden & Goossens 1991b;van der Linden et al. 1992), and cylindrical geometry (van der Linden & Goossens 1991a; Soler et al. 2011) under solar coronal conditions.Linear and follow-up nonlinear theory of thermal instability is a powerful tool to explain various fascinating features of the solar atmosphere.For example, the possible formation of a prominence in a current sheet is discussed by Smith & Priest (1977) and the dynamic thermal balance in a coronal arcade is studied in Priest & Smith (1979).The post-flare loop formation in a line-tied current sheet configuration with radiative energy loss was simulated in a 2D MHD setup by Forbes & Malherbe (1991).
More recently, multidimensional simulations related to prominence formation emerged in a variety of magnetic topologies.Xia et al. (2012) reported the ab initio formation of a solar prominence in a 2.5D MHD simulation in a bipolar magnetic arcade due to chromospheric evaporation and thermal instability.This was revisited in a quadrupolar arcade setup, in which reconnection was induced by the condensing prominence by Keppens & Xia (2014).That prominences can also form by feeding chromospheric matter within plasmoids during a flux rope eruption is developed by Zhao & Keppens (2022).3D models of prominence formation that establish the needed plasma cycle between chromosphere and corona are shown in Xia & Keppens (2016).More recently, prominence formation due to levitationcondensation (Kaneko & Yokoyama 2015;Jenkins & Keppens 2021) was demonstrated where a 3D realization is needed to allow magnetic Rayleigh-Taylor instability (Jenkins & Keppens 2022).The effect of thermal instability has also been explored for the formation and dynamics of coronal rain in magnetic arcades in 2.5D geometry (Fang et al. 2013(Fang et al. , 2015a)), in a weak magnetic bipole in 3D geometry (Xia et al. 2017), in a more self-consistent 3D radiative-magnetohydrodynamic setup in Kohutova et al. (2020), and for randomly heated arcades by Li et al. (2022) in a 2.5D geometry.
These works also triggered a renewed interest in more idealized studies of linear thermal instability and its nonlinear evolution, and in how the various linear MHD waves and instabilities may interact.Numerical analysis in the linear and non-linear domains due to the interaction of the slow MHD and entropy (thermal) modes was carried out in recent studies by Claes & Keppens (2019) and Claes et al. (2020), while the effect of different radiative loss functions on the onset and far nonlinear behavior of thermal modes was analyzed by Hermans & Keppens (2021).However, the influence of thermal instability on the tearing mode of solar current sheets has not gained much attention to date.Linear analysis by Ledentsov (2021a,b,c) shows that the instability growth rate in a pre-flare current sheet is modified if the nonadiabatic effects of radiative energy loss, resistivity and thermal conductivities are included.Sen & Keppens (2022)[SK22 hereafter] extended this into the non-linear domain and incorporated background heating and optically thin radiative loss into a series of 2D resistive MHD simulations.This study finds that the instability growth rate of tearing modes in a solar current sheet increases by an order of magnitude when these non-adiabatic effects are incorporated, such that we can meaningfully speak of coupled tearing-thermal evolutions.The 2D current sheet produced a chain of plasmoid-trapped condensations with cool material, which are thermodynamically similar to prominence (or coronal rain) in the solar atmosphere.
In this work, we extend our study to a 3D geometry and explore the tearing-thermal evolutionary process of an idealized current sheet model in solar atmosphere which is essentially non-adiabatic (with background heating, optically-thin radiative loss, and thermal conduction).Due to the mutual reinforcement of these instabilities, we demonstrate the combined effect of the complex evolution of the current layer, which disintegrates into finer structures with subsequent development of flux ropes, along with the formation of cool plasma condensations in the vicinity of this evolving current sheet.These localized cool, and plasma-condensed regions share similarities with the prominence and coronal rain structures observed in the solar atmosphere.Our findings here augment the growing theoretical basis for the combined effect of a current sheet fragmentation, and formation of cool-condensed plasma due to coupled tearing-thermal instability.This multi-mode evolution of the system occurring in association with or during reconnection in a current sheet is an important aspect to understand the dynamics and multi-thermal processes in the solar atmosphere.
The paper is organized as follows.In Sect.2, we describe the numerical model with an initial configuration with a precise magnetic and thermodynamic structure, and detail algorithmic aspects and boundary conditions.In Sect.3, we discuss the main results of the study, and relevance of the model in the solar atmosphere.Section 4 addresses the significance and novelty of the work for a typical coronal atmosphere, the scope for further development and points out how this work will be useful for future studies.

Numerical setup
We construct a resistive MHD simulation using MPI-parallelised Adaptive Mesh Refinement Versatile Advection Code or MPI-AMRVAC1 (Keppens et al. 2012;Porth et al. 2014;Xia et al. 2018;Keppens et al. 2021;Keppens, R. et al. 2023) in 3D Cartesian geometry.The spatial domain of the simulation box spans from −10 to 10 (in dimensionless units) along x − y − z directions.We activate the adaptive mesh refinement (AMR) up to level three, which gives a maximum resolution of 512 3 .If the box size is set in units of 10 Mm, this achieves the smallest cell size of 390 km in each direction.Automated refinement and derefinement is triggered based on the errors estimated by the instantaneous density (gradient) at each time step.
To explore the influence of thermal instability on the tearing mode in a current sheet, the following normalized MHD equations are solved numerically, Note that we use magnetic units where the magnetic permeability is unity.Here, I is the unit tensor, and ρ, T, B, v, and η represent mass density, temperature, magnetic field vector, velocity vector, and resistivity, respectively.A uniform resistivity, η = 0.001 (or 1.2 × 10 14 cm2 s −1 in physical units) is taken throughout the entire simulation domain.We adopt the Spitzertype thermal conductivity, κ || = 10 −6 T 5/2 erg cm −1 s −1 K −1 , which is purely aligned along the magnetic field.The total pressure p tot is the sum of the plasma and magnetic pressure given by where p is the gas pressure linked with the thermodynamic quantities through the ideal gas law.The total energy density is where γ = 5/3 is the ratio of specific heats for a monoatomic gas (fully ionized hydrogen plasma).We set up a current sheet configuration using the magnetic field components where B 0 = 1 (corresponding to 2 G in physical units) is the magnetic field strength, which is comparable with the observations in the solar corona, where the field strength at the height of 1.05-1.35solar radius are reported between 1-4 G (Lin et al. 2004;Kumari et al. 2019;Yang et al. 2020).The unit plasma density, temperature, and length scales are set as ρ = 2.34×10 −15 g cm −3 , T = 10 6 K, and L = 10 9 cm, which are relevant for the solar corona.The initial width of the current sheet is set to l s = 0.5 (5 Mm in physical unit), which is comparable with the observed flare current sheet thickness (Li et al. 2018;Savage et al. 2010).The magnetic field configuration given by Eqns.(9-11) represents a force-free field, and the polarity reversal of the magnetic field occurs around the z = 0 plane, where the current sheet is.In line with a fully force-balanced equilibrium of the system, we use isothermal and isobaric conditions as the initial setup of the model.
The third term at the right-hand side (RHS) of Eq. (3) represents the radiative cooling in the optically thin corona, where Λ(T ) is the cooling function developed by Colgan et al. (2008) and extended to lower temperatures following Dalgarno & Mc-Cray (1972).The precise temperature dependence of Λ(T ) was shown in Figure 1 of our previous 2D simulation (SK22).To maintain the initial thermal balance between optically thin radiative loss and background heating, H bgr of the system, we prescribe a uniform, time-independent value, The motivation of using the above form is that the radiative cooling term at the initial state exactly compensates the background heating term, and the heating/cooling (mis)balance in the system occurs after a long term evolution triggered by the external perturbations (which is magnetic field in this study).Note that, this heating model is similar, although uniform in space unlike to our earlier study in SK22.However, the role of different heating models based on the power-law behaviour of magnetic field strength, and density on the thermal runaway, and condensation processes have been reported by Brughmans et al. (2022), which finds that the different heating models can change the evolution and morphology of the condensations.Therefore, how the different heating rates can change the thermal balance of our model can be interesting to study in future.With a homogeneous density ρ 0 = 0.2 (4.68 × 10 −16 g cm −3 in physical units) and isothermal atmosphere T 0 = 0.5 (0.5 MK in physical units) as initial condition, we have an initially uniform plasma-β = 0.2, less than unity as appropriate for solar corona.We use the initial temperature, T 0 = 0.5 MK, which is in the temperature regime where the cooling function we use in this study has a very sharp gradient, and the heating-cooling mis-balance may be dominated due to perturbation from the equilibrium temperature.However, we also notice that the system achieves to thermal runaway state, and cool-condensed materials form for other equilibrium temperature regime, T 0 = 1 MK, which is shown in Appendix A. It is to be noted from the last term in the RHS of Eq. 3, that there is no role for thermal conduction at the initial time, as the system starts off isothermal.Therefore, the system is initially in thermal equilibrium, but the finite value of resistivity will drive the ideal force-balanced state away from its initial state, but only on the (slow) resistive timescale.This setup is liable to both linear resistive tearing modes, for which finite resistivity is key, and has all thermodynamic ingredients to allow for thermal instability.The equilibrium system is perturbed to trigger tearing modes, which can further trigger the thermal modes, and thus enforce each other in a coupled tearing-thermal fashion.We use parametrically controlled, monopole-free magnetic field perturbations mainly confined in the vicinity of the z = 0 plane (where the initial current sheet is present) and exponentially decaying for |z| > 0, Here, the parameter l = 20 × L matches the geometric sizes of the simulation domain along x and y directions respectively, the perturbation amplitudes ψ 01 = ψ 02 = 0.1 ensure a variation of 10% of the magnetic field strength B 0 , and we take the multimode distribution of the perturbations using n 1 = 4 and n 2 = 2.The magnetic field distribution (Eqs.9-11), and the perturbations (Eqs.13-14) ensure the solenoidal condition, After the initial setup, the system is allowed to evolve as governed by the Eqs.(1-6).The equations are solved numerically using a three-step Runge-Kutta time integration with a second order slope limited reconstruction method (Ruuth 2006) with 'Vanleer' flux limiter (van Leer 1974), and a Total Variation Diminishing Lax-Friedrichs (TVDLF) flux scheme.We follow the evolution of the system for up to 214.7 minutes and save the data at a cadence of 85.87 s, which gives 151 snapshots.We use periodic boundary conditions along x, y directions, and open boundary condition along the z direction.The wall clock time for the entire simulation run is ≈ 90 hours using 8 nodes with 288 processors in total with GNU-Fortran (version 6.4.0) compiler and open MPI 2.1.2.
Note that, there are some important differences in the initial conditions of this model with respect to SK22 mentioned in the following.(i) This is a force-free magnetic field configuration, and therefore an isobaric and isothermal medium ensures a fully force-balanced equilibrium, whereas the magnetic field configuration in SK22 was non-force free, and therefore we used a non-uniform density profile (though the initial temperature was also uniform) in such a way that it maintained the force-balance equilibrium, (ii) The magnetic field strength near the current sheet in SK22 is ≪ 1, which sets the plasma-β ≫ 1 near the current sheet, on the other hand, our current model has initially a uniform and low plasma-β = 0.2 throughout the entire simulation domain, (iii) the direction of imposed magnetic field perturbations for SK22 were both parallel and perpendicular to the current sheet, but in the current work we set the perturbation only parallel and concentrated around the current sheet plane.Besides the difference in affordable numerical resolution (SK22 has the maximum resolution of 2048 2 , whereas the current setup has maximum resolution of 512 3 ), the system studied here is intrinsically 3D, and is hence more relevant for actual current sheet conditions.

Global evolution
The spatial distribution of current density squared J 2 , as well as representative corresponding field line evolutions are shown in Fig. 1.The equilibrium configuration of the current sheet (at t = 0) is formed due to the fact that the magnetic field shears across z = 0 as given by Eqs. 9, 10, and 11.The initial current distribution is mostly oriented in the x−y plane and concentrated around z = 0, with its main contribution from J x = −dB y /dz and J y = dB x /dz (while J z is purely from the perturbed field).Due to the magnetic field perturbations, given by Eqs. 13 and 14, which are (mainly) confined near the z = 0 plane and extend along the x − y directions, linear resistive tearing modes start to develop leading to the disintegration of the current sheet.The inhomogeneities of J 2 in the current sheet plane due to the multimode perturbation (n 1 = 4, n 2 = 2) appear at the initial stage of the evolution as shown in Fig. 1a.Magnetic reconnections lead to pronounced magnetic topology changes, modifying the current sheet as shown in Fig. 1b.We see the perturbed field lines near z = 0 (see Fig. 1c) turn into extended flux rope-like structures while the system evolves (see Fig. 1d), yet the planar (x-oriented) field lines away from the current sheet plane remain unperturbed.The signature of flux ropes in the current density distribution in Fig. 1b can be clearly noticed.In a 2D setup, these would be the familiar magnetic islands or plasmoids, due to development of tearing modes in the current sheet.We notice the self consistent development of B z due to the perturbed fields around the current sheet plane as shown in Fig. 2a (note that the initial B z was set to zero with no magnetic field perturbation along z).Fig. 2b represents the plasma pressure distribution at the x − z plane (at y = 0), while the variation along the vertical green dashed line is shown in Fig. 2c, which shows the expected linear eigenmode structure of the tearing mode, seen as the kinks in the pressure distribution around z = 0.This implies the development of the tearing modes at the initial stage of the current sheet evolution.The nonlinearly developing tearing evolution creates density perturbations in the surroundings of the current sheet, and the radiative cooling (in combination with thermal conduction) becomes dominant over the constant background heating in localized regions of the domain.This triggers the cooling of those regions, which in turn condenses the regions even more.Hence, a runaway process starts where tearing and thermal modes reinforce each other, which causes the spontaneous growth of density and temperature inhomogeneities in the surroundings of the current sheet.Note that our box is periodic along x and y, so field aligned thermal conduction does not play a big role in the heating-cooling misbalance of the system (in the sense that it can not lead to heat fluxes down into lower-lying chromospheric regions: we have no stratification here).It just tries to homogenize the temperature along field lines, in competition with local resistive heating.
The temporal variation of the instantaneous maximal plasma density and minimal temperature are shown in Fig. 3a.The peak density increases from the initial uniform 4.68 × 10 −16 g cm −3 at t = 0 to 7.11 × 10 −14 g cm −3 at t = 214.7 min.A sharp rise of the peak density starts from ≈ 150 min.The instantaneous minimum temperature starts from 0.5 MK (the initial uniform equilibrium temperature), which also drops sharply at the same time where the density peak has the sharp rise (t ≈ 150 min), and reaches down to 10, 148 K at t = 214.7 min.This signals the formation of cool plasma condensations in the system at ≈ 150 min, or more than two-and-a-half hours following the initial tearing onset.The variation of the instantaneous maximal velocities (v x , v y and v z ) in Fig. 3b demonstrates that a dynamical instability of the system occurs at the condensation onset time (≈ 150 min).Here, the velocities are scaled in terms of Alfvén velocity, v a = 261 km s −1 , which is calculated based on the initial equilibrium density and the magnetic field strength of the system.We notice from Figure 4 that B z field evolves up to ≈ 0.4 G along the x − z plane (at y = 0) at t = 207.5 min, which is around an order of magnitude higher than the B z (x, z) field at t = 14.3 min shown in Figure 2a.The kinks appear in the bottom panel of Figure 4 in the B z field around z = 0 plane (at y = 0) along x = ±2 Mm shows the evolution of tearing mode around the current sheet plane.The evolution of density and temperature distributions along three orthogonal planes at x = 0, y = 0, and z = 0, and its 3D visualization are shown in Fig. 5. Density and temperature inhomogeneities appear in and around the current sheet plane that are aware of the initial multimode (n 1 = 4, n 2 = 2) magnetic field perturbation, as shown at t = 14.3 mins in Figs.5a and 5d, respectively.Due to tearing-associated thermal instability, localized condensed structures can be clearly seen in the x = 0, y = 0, and z = 0 plane at t = 207.5 min (see Fig. 5b).The 3D visualization of the density distribution is shown in Fig. 5c, where we see the condensed structures are formed around the current sheet plane.These condensed structures correspond to the cooler (∼ 10 4 K) regions compared to the background medium (see Figs. 5e and 5f).Their order 100 density and temperature contrasts are similar to coronal rain or prominence features.Note that they happen near the evolving current sheet, which is heated up to several million degrees due to effective Ohmic heating.The periodic boundary treatments that we use in this study are acceptable for the entire evolution, since the plasmoid sizes do not reach the lateral domain size of the simulation box.The condensations from thermal instability develop locally, and are expected not to be influenced by the type of boundary used laterally.Histograms of the mass and temperature distributions in the entire domain at t = 207.5 min are shown in Fig. 6, where we see in Fig. 6a that 98.8% of the cells within the simulation box contain mass within the range of 4.84 × 10 4 to 1.43 × 10 5 kg (the total number of cells in the simulation box used here is the effective resolution 512 3 ), while the number of cells with mass ≳ 1.43 × 10 5 kg is very low (≈ 1.2%).The mass range with the most number of cells is in accord with the mass determined by the initial equilibrium density of the medium (since density remains almost unperturbed away from the current sheet plane).Similarly, most cells contain the temperature value in the range of ≈ 0.32 − 0.63 MK, namely 75.1% of the total cells in the box as shown in Fig. 6b.Note that the initial equilibrium 0.5 MK temperature of the system, lies within this range.The fraction of the cell numbers with temperature ≲ 10 5 K, and ≳ 1 MK are low, and they are 2.1% and 1.2% respectively.This implies that the cool-condensed structures, as well as the hot regions with ≳ MK temperatures (which appear around the current sheet plane due to reconnection-induced heating) are very localised in the medium.To appreciate the thermodynamics of the cool condensed structures, we show slices of plasma pressure and different velocity components in Fig. 7.We notice that the velocities that develop in the different cutting planes in Figs.7d-7f are associated with pressure gradient driven flows (see Figs. 7a-7c), also called siphon flows.They are directed from higher to lower pressure regions and demonstrate sub-Afvénic speeds (Alfvénic Mach number reaches up to ≈ 0.075).Due to plasma accumulation as evident from the velocity maps, the cool condensation sites develop in these same regions as shown in Fig. 5.When the thermal runaway sets in after a long term evolutionary process, the velocities are domi-Article number, page 4 of 12 sites, and drag the magnetic field lines along, which entangle and form flux ropes.The heating/cooling (mis)balance in different planes are shown in Fig. 8.We use a uniform (and constant in time) background heating of 6.244×10 −53 erg g 2 cm −3 s −1 which is equal to the radiative loss (no field aligned thermal conduction at the initial state due to isothermal condition) at the initial equilibrium state of the system.But this balance breaks down due to tearing influenced thermal changes around the current sheet, and the radiative loss in some regions near the current sheet plane dominates over the heating.These regions correspond to the cool condensed structures as shown in Fig. 5.The regions away from the current sheet plane maintain the heating/cooling balance as the initial perturbation was concentrated around the current sheet plane, and therefore those regions maintain the initial equilibrium density and temperature of the system.The energetic evolution of the system is shown in Fig. 9.As expected (despite the open boundaries along z), this shows that the mean total energy density (E T ), which is a sum of mean kinetic (E k ), magnetic (E M ), and internal (E int ) energy densities (see Eqn. 8) given by respectively (where, V is the total volume of the simulation box) is nearly conserved in time.The resistivity and thermal conduction effects do not cause any deviation from total energy conservation, and only the heating/cooling misbalance may lead to net energy losses (or gains).Due to the resistive MHD evolution of the system, the energy exchange between magnetic and internal energies show an anti-correlation nature while the system evolves.This energy exchange occurs due to the mean Ohmic dissipation, which is quantified as We notice that the mean magnetic energy density decreases up to ≈ 150 min due to release of magnetic energy through reconnection.Thereafter, when the thermal runaway process happens in the system in a coupled tearing-thermal fashion, the field lines start to entangle with each other and form flux ropes.This generates magnetic stress, and leads to the enhancement of the magnetic energy density.However, the open boundaries along z−direction, which are sufficiently away from the central current sheet plane allow the magnetic energy flux to flow through the boundaries.The current density (spatially) distributes rapidly due to the implemented perturbation, and therefore we see a sharp drop of volume averaged J 2 initially.The temporal evolution of the mean kinetic energy density shows that the instability dynamics of the system shows a (nearly) quasi-equilibrium nature up to ≈ 150 min (which is the onset time of condensations), and then rises sharply due to a tearing-thermal coupled unstable evolution.

Synthetic Observation
The Solar Dynamics Observatory Pesnell et al. ( 2012) is a near-Earth orbiting satellite suite capable of routinely observing the Sun from its photosphere to corona.The Atmospheric Imaging Assembly Lemen et al. (2012) on board captures images of the solar atmosphere with a temporal cadence of ∼ 12 s at a spa-tial resolution of ∼ 0.6 ′′ per pixel across a range of ultraviolet and extreme ultraviolet wavelengths, the latter mainly associated with the different ionisation states of Iron, namely Fe xiixxiv.Emission from such highly ionised Iron corresponds to coronal temperatures in the broad range of a few hundred kK to around 20 MK.The emission coefficient for such coronal plasma is given by, where A b is the abundance of the emitting species, n e is the ambient electron number density, and G λ is the contribution function for a specific wavelength, indicated to be additionally dependent on n e as well as temperature T .In the absence of a modelled n e , it is instead approximated using local thermodynamic equilibrium Saha-Boltzman.This contribution function is precomputed for each of the Atmospheric Imaging Assembly passbands using the CHIANTI atomic package for a range of electron number densities and temperatures between 10 6 -10 12 cm −3 and 10 4 -10 8 K, respectively (Landi & Reale 2013;Verner et al. 1996).τ here denotes the local optical depth, computed within each voxel that a given line of sight intersects as the product of the local absorption coefficient α λ and the length of the ray within that voxel.
The atmosphere of the Sun is optically thin at the wavelengths corresponding to the emission by these Fe lines.As such, the standard approach to synthesizing simulations so as to resemble the appearance of structures as seen by Solar Dynamics Observatory/Atmospheric Imaging Assembly is to employ Eq. 19 in a local manner and apply an arbitrary line of sight integration according to the position of an observer.For structures that are majority-comprised of material at coronal temperatures, such as coronal loops, this is deemed a sufficient approach (Van Doorsselaere et al. 2016; Gibson et al. 2016).
Cool condensations within the solar corona appear dim in extreme ultraviolet contrast (cf.Carlyle et al. 2014), and indeed the value of G λ for the Atmospheric Imaging Assembly passbands that we consider is many orders of magnitude lower at condensation temperatures than for coronal temperatures.However, the strong contrast is not only due to small G λ , but also the direct ab-Fig.9: Time series of mean kinetic (E k ), magnetic (E M ), internal (E int ), total (E T ) energy density, and ohmic heating rate (E ohm ), normalised with respect to their maximum values which are 1.89 × 10 −3 , 1.61 × 10 −1 , 7.40 × 10 −2 , 2.13 × 10 −1 erg cm −3 , and 2.23 × 10 −9 erg cm −3 s −1 , respectively.sorption and removal of background EUV photons from the light beam (Kucera et al. 1998).This is due to a number of the wavelengths observed by Solar Dynamics Observatory/Atmospheric Imaging Assembly lying below the head of the Hydrogen Lyman continuum at 912 Å, and so H, He, and He i (with characteristic temperatures < 10 kK) are photo-ionised by this extreme ultraviolet emission up to the ionisation continuum of He ii at 227 Å (Williams et al. 2013).Hence, extreme ultraviolet photons are progressively removed from the line of sight if such cool material is encountered.The absorption coefficient as a consequence of this extreme ultraviolet photo-ionisation can be approximated in local thermodynamic equilibrium by, where s refers to the photo-ionised element and A b,s and σ s are the assumed abundance and cross-section of ionization of element s, measured observationally and experimentally/theoretically, respectively.The summation weights w s (τ) are the ratios of the number densities as w H = 1 − n Hi /n H , w He = 1 − (n Hei − n Heii )/n He , and w Hei = n Heii /n He .To obtain an approximation to these weights, and in accordance with our previous assumptions of local thermodynamic equilibrium, we iteratively solve for each of the considered population densities using the Saha equation and associated partition functions.We consider convergence under the local thermodynamic equilibrium assumption to be achieved once the absolute relative difference of n e between iterations drops below an arbitrary value of 10 −4 (a method developed by Zhou et al. 2019).One must necessarily consider this photo-ionisation to correctly approximate the appearance of cold plasma condensations, if present, when synthesizing simulations of the solar atmosphere (Jenkins & Keppens 2022).
Both the emission and absorption quantities as defined above are purely local properties, the total emergent intensity I λ (τ λ ) along a given line of sight through these local voxels is then given by the integral form of the transport equation, where the combined influence of j λ and α λ are taken into account in the source function S λ = j λ /α λ , and τ λ now represents the total optical thickness along the chosen line of sight, and I λ (0) is the intensity of any background illumination, when the emergent intensity is calculated along the specific LOS with zero optical thickness region (Rybicki & Lightman 1986).The non-standard inclusion of the absorption coefficient requires every local voxel, with their respective optical depths, τ ′ λ to have access to a globally integrated and line of sight-specific optical depth, τ λ .Such a requirement is not compatible with the block-based architecture of MPI-AMRVAC and is hence completed in post processing using a combination of yt-project, numpy, scipy, and matplotlib in python.The implementation here represents an update of that previously presented in Jenkins & Keppens (2022).
Ground-based observatories do not have access to extreme ultraviolet wavelengths as recorded by Solar Dynamics Observatory/Atmospheric Imaging Assembly, and instead commonly observe, amongst others, the strong Hydrogen n = 3 → n = 2 (Hα) line at 6563 Å.This line is known to straddle the opticallythick -thin divide under solar atmospheric conditions, and a complete handling of the plasma-light interaction of such photons through the simulation domain would require non-"local thermodynamic equilibrium" modelling, outside the scope of this study (but we point an interested reader to the recent work of Jenkins et al. 2023, for comparison).Instead, Heinzel et al. (2015) reported an approximate approach to relating local pressure and temperature conditions to Hα opacity (α λ ) according to their series of 1.5D radiative transfer models.A key property of these models considers the source function of Eq.21 to remain constant along a given line of sight, and enables the following simplification, The resulting emergent (specific) intensity of the Hα line is therefore found through a line of sight integration of the approximate α λ according to the tables of Heinzel et al. (2015), wherein the authors also provide a coarse height-dependent estimation to the constant value of S λ .
To convert the physical variables (plasma density and temperature) into spectroscopic observables (namely specific intensity), we generate the synthetic maps of the simulation output using forward modeling.Measuring LOS integrated specific intensity depends on the (theoretical) viewing position of the observer.Thus, due to the spatial distribution of the substructures due to the condensations, the synthetic maps show differences for different LOS views.Fig. 10 shows the synthetic maps of the reoriented spatial domain (we now show the current sheet vertically) capturing the internal structures as being viewed sideways, representative for the solar limb.We do so for three different extreme ultraviolet passband filters of AIA, 171, 193, and 304 Å using the contribution functions of their specific spectral lines, which each highlight material at ≈ 0.8 MK, 1.5 MK, and 80 kK respectively, for a LOS integration along our y direction at t = 207.5 min.An animation of the synthetic maps for different LOS directions around the z axis is available online.Due to the absorption features of the condensed plasma regions, the strongest absorption corresponds to the LOS direction along the y direction, as the condensations are aligned with it (see Fig. 5).The cool, condensed plasma appears darker in AIA 171, 193, and 304 Å passband filters due to photo-ionisation of H, He, and He i as previously outlined.For the remaining optical Hα line, we obtain a positive intensity in the absence of any background illumination, as would be the case for, amongst others, prominences and jets positioned at or above the limb.The island like structures near z = 0 in the EUV maps in Fig. 10   The resolution of the presented simulation is 390 km 2 per pixel.In all cases, the spatial resolution of our synthesized images have been modified to match that of an equivalent observatory.For the EUV pass bands, this is set to the instrumental resolution of the AIA filters (≈ 430 km 2 per pixel), and leads to a minor smearing of the finer structures inside the finest condensations.This blurring is stronger for the Hα line where the resolution is set to 1 ′′ = 725.27km so as to match the GONG ground-based instruments (Harvey et al. 2011).It will be useful to increase our model resolution for comparisons against the state-of-the-art observations anticipated from Solar Orbiter High Resolution Imager (HRI) campaigns (Rochus et al. 2020).

Summary and Conclusions
The study of the combined tearing-thermal instability of an idealized 3D current sheet configuration as addressed in this work is to understand the theoretical basis of the multi-mode evolution of a current sheet, which is an important aspect to understand the dynamics and multi-thermal behaviour of solar atmosphere.In contrast to our earlier 2D simulations of a non-force-free current sheet (SK22), where thermal runaway was happening simultaneously with chaotic tearing, and condensations were trapped and gathered within coalescing plasmoid structures, the current setup shows a clear tearing evolution at first leading to 3D topological changes in the magnetic field, and later on show runaway condensations near the central current sheet.Note that there are some important differences between the initial setup of the current model with the one used in SK22, as explained in section 2, most notably perhaps the plasma beta regime which is uniformly low in the current 3D simulation.Here, the condensation onset time in the current model is much later than in SK22, because the initial setup and adopted perturbations do not directly modify the heat-loss balance and a purely tearing evolution starts at first.In a 3D long-term simulation of a macroscopic current sheet, we find that cool plasma condensations are produced in the vicinity of the current sheet due to tearing-influenced thermal instability.Our findings are based on a 3D resistive MHD simulation with non-adiabatic effects of radiative cooling (for an optically thin medium), background heating, and magnetic field-aligned thermal conduction, which are relevant for the solar corona.We find that the plasma density of the condensations can go up to ∼ 10 −14 g cm −3 , which is two orders of magnitude more than the initial background density (4.6 × 10 −16 g cm −3 ), and the temperature of the condensations can drop down to ∼ 10 4 K, which is an order of magnitude less than the initial equilibrium temperature (0.5 MK) of the medium.These locally, in-situ forming cool condensed structures hence show similar thermodynamic contrasts with their surroundings, just as coronal rain or prominences observed in the solar atmosphere.This is highlighted by synthetic views, which account for important absorption effects within the synthesis of the EUV channels of Solar Dynamics Observatory/Atmospheric Imaging Assembly.However, our model ignores the effect of the stratified solar atmosphere due to solar gravity.The condensation time scales in this idealized current sheet model is ≈ 150 min, which is larger than the time scales (which is ≈ 30 min) compared to the earlier study for a post-flare coronal rain model by Ruan et al. (2021), where they use a stratified atmosphere, and reveal the multi-thermal aspects of a post flare loop underneath the current sheet and reconnection sites.Whereas, our study demonstrates the development of multi-thermal plasma (∼ 10 kK -MK) in and around the current sheet and reconnection sites.Also, earlier models of coronal rain in arcades show a clear tendency to develop strong shearing motions (Li et al. 2022;Zhou et al. 2021;Fang et al. 2015b), and velocity shear can alter the tearing mode growth, which we plan to explore by incorporating velocity shear flows into the model.Future efforts should exploit a more realistic 3D model by using the initial magnetic field configuration from an extrapolated magnetic field for an active region vector magnetogram.Nevertheless, the current study sheds new light on the instability of solar current sheets due to the combined effect of tearing and thermal modes, and unifies multi-thermal processes in current sheets with the formation mechanism of cool condensations such as prominences, and coronal rain in solar atmosphere, which are important aspects in the understanding of the broader solar coronal heating.

Fig. 1 :Fig. 3 :
Fig. 1: Isosurface (five isosurface values are taken from instantaneous minimum to maximum) plots of the current density squared J 2 , which represent the 3D structure of the evolving current sheet.Panel (a) represents the perturbations appearing in the current sheet plane due to the multimode (n 1 = 4, n 2 = 2) magnetic field perturbations in Eqns.(13-14), (b) represents a time when the current sheet disintegrated due to tearing-thermal coupled evolution of the system.Panels (c) and (d) represent the magnetic field corresponding to panels (a) and (b), respectively.Distances along x, y, and z axes are scaled in units of 10 4 km, where we see the field lines perturbed at 14.3 min around the current sheet plane only, which evolves to sizeable flux rope structures at 207.5 min.(An animation of the figures is available online.)

Fig. 5 :Fig. 6 :
Fig. 5: Spatial distribution of density (top panel) and temperature (bottom) in the 3D domain, where the distances along x, y, and z directions are in units of 10 4 km.Density (a and b), and temperature (d and e) distribution along three orthogonal slices along x = 0, y = 0, and z = 0 planes for two different times, t = 14.3 and 207.5 min are shown.(c) and (f) represent isosurface views (five isosurfaces ranging from minimum to maximum values) on density and temperature, respectively.Panels (a) and (d) are the early phase of the evolution, where the density and temperature inhomogeneities appear around the current sheet plane due to the multi mode magnetic field perturbation.Panel (b) and (c) illustrate where high density structures appear, cospatial with cool (∼ 10 4 K) regions in (e) and (f).(An animation is available online.)

Fig. 7 :
Fig. 7: Distribution of plasma pressure (top panels) and velocity (bottom panels) at t = 207.5 min in different planes, where the velocities scaled in units of Alfvén velocity, v a ≈ 261 km s −1 .
represent plasmoids, which are the manifestation of the extended flux ropes

Fig. 10 :
Fig. 10: Specific intensity counterparts to the simulation, where we account for emission and absorption.Synthetic maps for the broadband 171, 193, 304 Å SDO/AIA filters, and a narrowband Hydrogen-Hα filter, for a LOS view along the y direction at t = 207.5 min, are shown from the left to right panels, respectively.An animation of this figure for a rotating LOS view around the z (from 0 to 360 • ) axis is available online.This shows the limb view on a current sheet that shows thermal-tearing evolutions.