Two-fluid reconnection jets in a gravitationally stratified atmosphere

The density decreases exponentially with height in the solar gravitationally stratified atmosphere, therefore the collisional coupling between the ionized plasma and the neutrals also decreases. Here, we investigate the role of collisions between ions and neutrals on the reconnection process occurring at various heights in the atmosphere. We perform simulations of magnetic reconnection induced by a localized resistivity in a gravitationally stratified atmosphere, where we vary the height of the initial reconnection X-point. We compare a magnetohydrodynamic (MHD) model and two two-fluid configurations: one where the collisional coupling is calculated from local plasma parameters and another where the coupling is decreased, so that collisional effects are enhanced. Simulations in a stratified atmosphere show similar structures in MHD and two-fluid simulations with strong coupling. However, when collisional effects are increased to attain representative parameter regimes, we find a nonlinear runaway instability, which separates the plasma-neutral densities across the current sheet (CS). With increased collisional effects, the initial decoupling in velocity heats the neutrals and this sets up a nonlinear feedback where neutrals migrate outside the CS, replacing charged particles which accumulate towards the center of the CS. The reconnection rate has a maximum value around 0.1, similar for both reconnection heights, and is consistent with the use of a localized enhanced resistivity used in all three models. The initial stages of plasmoid formation, observed near the end of our simulations, is influenced by the outflow from the primary reconnection point, rather than by collisions. We synthesize optically thin emission for both MHD and two-fluid models, which can show a very different evolution when the charged particle density is used instead of the total density.


Introduction
The temperature increases steeply towards the corona of the quiet sun, throughout the transition region located at ≈ 2.5 Mm, and coronal plasma becomes fully ionized (O'Flannagain et al. 2015).Observations showed that field lines below ≈ 2.5 Mm can change their connectivity in about 1.5 h, suggesting fast reconnection mechanisms in the solar chromosphere (Close et al. 2004;Yamada et al. 2010).
The reconnection mechanism involved in solar flares releases a huge amount of magnetic energy, therefore it is believed that the reconnection is driven at large scales, evolving towards a Petschek-type configuration (Sato & Hayashi 1979).Petschek (Petschek 1963) proposed a stationary 2D reconnection model that extends a Sweet-Parker (Sweet 1958;Parker 1957) layer (the diffusion region) at the center with standing slow-mode shocks from its ends.The Sweet-Parker reconnection rate scales as ∝ S −1/2 , where the nondimensional Lundquist number, S is defined as S = v A L/η, with v A being the Alfvén speed, L a hydrodynamic characteristic length and η the value of the resistivity (in units of m 2 /s) This classical resistive MHD model proposed by Petschek (1963) predicts reconnection rates much higher than the Sweet-Parker model, where it scales as ∝ ln −1 (S).For a value of S = 10 6 , the Petschek-type reconnection gives a reconnection rate 65 times larger than the Sweet-Parker model (Leake et al. 2012).With larger resistivity, the slow shocks become wider (Sato & Hayashi 1979).The magnetic energy decreases, being converted into thermal energy through Ohmic heating near the diffusion region, and to kinetic energy through the work done by the Lorentz force, away from the diffusion region, especially at the location of the slow shocks (Ugai 1995).In simulations of driven reconnection, anomalous resistivity is assumed in order to impede the unbounded growth of the current density (Sato & Hayashi 1979;Ugai 1995;Uzdensky 2003).
Anomalous resistivity can be produced by ion acoustic turbulence, and is then a function of the ion electron drift velocity, as shown by Manheimer & Flynn (1971).A more simplified, but frequently adopted model of anomalous resistivity is a space localized resistivity around the X-point, which is also used in this paper.Unlike anomalous resistivity which depends on the current density, the localized resistivity only depends on space, but both prescriptions lead to fast magnetic reconnection (Biskamp & Schwarz 2001).This locally enhanced resistivity leads to fast reconnection, however, it is not very clear whether fast reconnection is due to a high value of the local resistivity in the diffusion region, or to the localization (Biskamp & Schwarz 2001).Numerical 2D resistive MHD simulations found that having a flat local resistivity profile near the X-point can induce spontaneous symmetry breaking in the otherwise symmetric Petschek configuration Baty et al. (2014).In this paper, we will have an additional broken up-down symmetry due to gravity from the beginning, and extend Petschek-reconnection findings to a two-fluid, plasma-neutral setting.
Many simulations study the standard solar flare scenario, where a vertical current sheet (CS) evolves to form post-flare loops.Takasao et al. (2015) studied post-flare loops in the MHD approximation using a 2.5D setup, where gravity is neglected and thermal conductivity is adopted along the field lines.In their simulations, Petschek-type reconnection develops because of localized resistivity and the slow shocks are essentially isothermal due to effective thermal conductivity.Wang et al. (2021) study reconnection in a 2.5D MHD setup, using a spatially localized resistivity and an analytic density profile, with uniform pressure.In their model, gravity is neglected, while anisotropic thermal conductivity is incorporated.The reconnection rate was found to be slightly smaller when plasma β increases, and the rate is also smaller with thermal conductivity.The reconnection rate reaches a maximal value of 0.01.The authors show that reconnection at the termination shock due to interaction between magnetic islands formed along the primary current sheet (CS) and the magnetic arcade below is almost as important as reconnection in the main CS for releasing magnetic energy.Jets produced by MHD simulations of Petschek reconnection in a 2D setup without gravity have properties of small-scale flares observed in the solar atmosphere (Innes & Tóth 1999).State-of-the-art solar flare simulations extend these efforts with the inclusion of gravitational stratification, and even include the effect of fast electron beams that self-consistently interact with large scale MHD simulations, identifying many ingredients found in actual observations (Ruan et al. 2020).The step to full 3D MHD simulations, including gravity, thermal conduction and reproducing turbulent regions consistent with observed non-thermal velocities was made in Ruan et al. (2023).EUV synthetic images produced from these 3D flare models show very good agreement with observations.However, extension to plasma-neutral setups are needed to address more chromospheric flare counterparts, and this work is a first step towards that goal.
Reconnection jets have been observed at all heights in the solar atmosphere, from photosphere to corona (Schmieder et al. 2022) in both cool and hot spectral lines.Takasao et al. (2013) suggest that spicules are chromospheric jets in emerging flux regions, which disappear in chromospheric line images before returning, probably due to heating.Hot and cool jets observed by Schmieder et al. (1995) were suggested to form in the lower corona or upper chromosphere.Aspects of coronal X-ray jets were successfully reproduced in simulations by Yokoyama & Shibata (1996).Chromospheric anemone jets with velocities of 10 km/s comparable to the local Alfvén speed were observed in the upper chromosphere, but could not be observed in the lower chromosphere, where the Alfvén speed is much lower (Takasao et al. 2013).
Many observations show clear signatures of plasmoids formed during the reconnection process, and track their motions using emission signatures.Plasmoids have been observed as periodic blobs in optically thin AIA lines (Kumar et al. 2023).These plasmoids can appear in the nonlinear evolution of current sheets that are liable to linear resistive tearing modes.Since plasmoids form on the current sheets that also develop Petschek-type configurations with outflows, once they are formed, it is relevant to study the linear stability of a CS in the presence of outflows.In an early 2D purely linear MHD analytical study, it has been shown that outflows have a stabilizing effect for the tearing mode (Bulanov et al. 1978).In this paper, we will discuss stability aspects of a CS due to tearing in a two-fluid setting.Our simulations contain plasmoids, and we will make synthetic images that directly relate to the observed blob features.
Since our work is using a plasma-neutral two-fluid model, several works that looked into two-fluid reconnection are of direct relevance to our study.2D simulations in the two-fluid approach show that the reconnection rate is increased because of recombination and larger outflows (Leake et al. 2012).Ionization/recombination processes would put additional constraints on the background stratification, leading to nontrivial equilibrium conditions (Snow & Hillier 2021).A non-static equilibrium introduces a new free parameter through the gradient in the vertical velocity, moreover it can explain the formation and the properties of the transition region (Song et al. 2023).In this paper, we generalize the work to stratified settings, but will not include ionization/recombination processes in our simulations.Leake et al. (2013) obtain high reconnection rates around 0.1 due to two-fluid effects, which otherwise would be obtained by using Hall, kinetic effects or localized resistivity.In stratified setups that were liable to the Rayleigh-Taylor instability, secondary reconnection events showed that two-fluid effects are locally very important (Popescu Braileanu et al. 2023).Murtas et al. (2021) found that plasmoid coalescence happens faster in the two-fluid model than in MHD, because the effective Alfvén speed based on a twofluid density defined by the collisional coupling, is larger.In 1D two-fluid simulations of slow magneto-acoustic shocks (of the type relevant in a Petschek-type reconnection region), frictional heating leads to a localized region around the 'reconnection' point with increased temperature in the neutral fluid Hillier et al. (2016).As a consequence a blast wave in the neutral fluid develops with overshoots in the neutral density and velocity (Hillier et al. 2016;Snow & Hillier 2019).We will show in our multi-dimensional setup how the detailed CS structure may show intricate decoupling (and runaway) effects in the Petschek-type configuration.Here, we extend previous works by studying reconnection in a 2D stratified setup, accounting for the presence of coupled plasma-neutral species.
We present the numerical setup in Section 2, the results of our simulations in Section 3. We then create synthetic views from the simulation snapshots, presented in Section 4 and we summarize our conclusions in Section 5.

Numerical setup
We consider a gravitationally stratified atmosphere where we define a temperature profile with height z described by: where w tr = 0.2 Mm, the transition region height is z tr = 2 Mm, the chromospheric and coronal temperature is set through T ch = 8 × 10 3 K and T co = 1.8 × 10 6 K.The initial   2. Collisional mean free paths as functions of height.We show the mean free paths between ions and neutrals (λin, red lines) and between neutrals and ions (λni, blue lines) for the two collisionality regimes considered: self-consistent from plasmaneutral parameters (2fl, solid lines, left axis) and reduced case 2flα (dashed and dotted lines, right axis, note the different order of magnitude).The collisional mean free paths are calculated according to Eq. ( 8).
profiles of temperature and the neutral and charged, as well as total fluid densities are shown in Figure 1.We use an ideal equation of state for both neutrals and charges and the normalized mean molecular weight is considered uniform and constant, having the values of 0.5 for charges and 1 for neutral species (Popescu Braileanu & Keppens 2022).
We use a force-free sheared magnetic field, changing its far-field (large horizontal coordinate |x|) vertical orientation near x = 0, with uniform magnitude B 0 = 10 −3 T , given by where L s = 0.02 Mm.A similar force-free profile has been used in the simulations of Ruan et al. (2020Ruan et al. ( , 2023)).The current sheet width calculated as the full width at half maximum (FWHM) of the corresponding current density component J y0 is 0.036 Mm.We will use a localized resistivity of the particular form where η 0 ≈ 8 Ωm and η 1 ≈ 0.8 Ωm.The reconnection point will be at x = 0 always, but we will vary the reconnection height z rec , as mentioned in Section 2.2 below.To trigger the reconnection, we adopt an initial velocity variation.The initial perturbation is in the x-direction only, trying to bring the field lines closer around the reconnection point, having the form: with and the total Alfvén speed, having ρ tot = ρ n + ρ c the total density.We choose the amplitude A = 10 −1 .In a twofluid simulation, this initial perturbation is the same for the velocity of charges and neutrals.
The two-fluid model uses the newly developed module in the fully open-source MPI-AMRVAC code (Popescu Braileanu & Keppens 2022; Keppens et al. 2023).The equations solved are the nonlinear, compressible, resistive two-fluid Eqs. ( 1)-( 7) from (Popescu Braileanu & Keppens 2022).In a 2.5D geometry the domain xz with x ∈ [−0.5, 0.5] Mm and z ∈ [0, 7] Mm is covered by a grid with base resolution of 256×1024 and we use five levels of refinement, having an effective resolution of 1024×4096 points and the size of the finest cell ∆x=9.76×10−4 Mm and ∆z=1.7×10−3 Mm.The bottom boundary in the z-direction is closed (antisymmetric for vertical velocities and symmetric for the rest of the variables) and we use open boundary conditions (symmetric for all the variables) for the top boundary and both side boundaries in the x-direction.The region −0.2 ≤ x ≤ 0.2 is always refined at the highest level, so that the CS is properly resolved.The refinement criterion is based on density only for the MHD cases and equally on charged and neutral density in the two-fluid runs.We use the splitting of the equilibrium force-free magnetic field (Xia et al. 2018) and the gravity stratification for both neutrals and charges (Yadav et al. 2022;Popescu Braileanu & Keppens 2022).In this approach, the magnetic field B, densities ρ n , ρ c and pressures p n and p c are split into a time-independent, B 0 , ρ n0 , ρ c0 , p n0 , p c0 and time-dependent B 1 , ρ n1 , ρ c1 , p n1 , p c1 parts.The equations solved in the code are for timedependent quantities, while the equilibrium conditions are explicitly removed from the equations: Mathematically, the split equations are equivalent to the unsplit equations, but numerically, the splitting helps avoiding an unwanted evolution due to numerical dissipation of the equilibrium.

Coupling aspects
Because of the very low mass of electrons compared to ions, the collisions between charges and neutrals are effectively collisions between ions and neutrals, and the mean free path between ions and neutrals and between neutrals and ions can be defined as: where ν in = αρ n and ν ni = αρ c denote collision frequencies between ions and neutrals and between neutrals and ions, respectively (see Eq. ( 15) in Popescu Braileanu & Keppens 2022).The characteristic velocity is the Alfvén speed of the whole fluid, as defined (generalized to v A (x, z, ; t)) by Eq. ( 6) (see also Eq. ( 37) in Popescu Braileanu & Keppens 2022, where the characteristic velocity was chosen as the fast magneto-acoustic speed).The collisional parameter α (see Eq. (A.3) in Popescu Braileanu & Keppens 2022) is defined as: where the collisional cross-section considered here is Σ in = 10 −19 m 2 .T cn is the average of temperatures of the neutrals and charges.In this paper, we consider two cases for the collisional coupling, one when α(x, z; t) is calculated consistently from instantaneous plasma parameters and another where it is set to a constant value.When α is calculated self-consistently, its profile varies slowly with height, with values initially between 5.8×10 12 m 3 /kg/s and 8.2 × 10 12 m 3 /kg/s.These minimum and maximum values remain practically unchanged at the end of the simulation.These high values imply that the coupling is near perfect (and hence MHD-like behavior is expected) throughout the domain.We will compare this to a run where we instead set α to a constant value throughout.This constant value of α = 3.84 × 10 8 m 3 /kg/s is almost four orders of magnitude smaller than the consistently computed values.However, as we argue in the section below, this reduced coupling value ends up to be more representative for actual solar settings.Moreover, in that regime, we will have two-fluid effects that are important in stratified reconnection setups.

Parameters of the simulations
For our study of reconnection in stratified settings, we will compare three collisional regimes, namely: (1) a single fluid MHD model (label "MHD"); (2) a two-fluid plasma-neutral model where the collisional parameter α is calculated selfconsistently from plasma values (label "2fl"); (3) a two-fluid model where the collisional parameter α is constant and set to a smaller value of 3.84 × 10 8 m 3 /kg/s (label "2flα").
We will study 6 cases in total, since each collisional regime will vary the reconnection point from z rec = 2 Mm, with reconnection in the upper chromosphere to low transition region, to z rec = 4.5 Mm, for coronal reconnection.
In the MHD model, the initial equilibrium atmosphere is constructed by summing the densities (shown in Figure 1) and pressures for charged and neutral fluid at all heights.The resulting mean free paths between ions and neutrals and between neutrals and ions for the two two-fluid models, 2fl and 2flα, are shown in Figure 2. In this 2flα case, both values of the mean free path in the transition region and above are O(0.01Mm),hence they are similar to the width of the CS, while in the 2fl case the values of the mean free path are O(10 −6 Mm).Because of the weak dependence on height of the value of α calculated self-consistently, both profiles of the mean free paths (2fl and 2flα) look similar, being dominated by the variation of the density, included in the calculation of the mean free path as from Eq. 8.The mean free path does not change significantly during the simulation, being rather determined by the equilibrium plasma parameters.
We will now argue that the 2flα case is more solarrelevant.Indeed, we used a temperature profile defined by Eq. 1, which is slightly smoother than VALC (Vernazza et al. 1981), but it has the advantage that we directly control the width (w tr ) of the transition region (TR).A similar temperature profile defined by an analytic function has been used by other authors (Dover et al. 2021).This overly smooth temperature variation also implies a larger fraction of neutrals at both used reconnection heights, especially at the coronal reconnection point z rec = 4.5 Mm, taking into account that the scale height of the neutrals is twice smaller than that of the charged particles.The other reconnection case has z rec = 2 Mm, i.e. starts its reconnection at the middle of the TR.How effective the collisional coupling between plasma-neutral species is at these heights, is largely set by the densities they attain there.
When we integrated the vertical profiles, at the base of the atmosphere at z = 0, we used the total number density from the VALC model namely n T ≈ 10 23 m −3 , but we had to consider an ionization fraction of 0.1, so that we still obtain more charges in the corona despite the smoothing of the temperature profile.Hence, the adopted temperature variation, together with the imposed bottom densities lead to a reversal of the dominance of neutrals over charges at z = 1 Mm, while the entire transition region and corona is charge-dominated.We find that the overly smoothed temperature profile actually increased the density in the upper chromosphere and corona above observed values, making collisional coupling now larger there.This justifies the use of a smaller and more representative value of α in the simulation 2flα, which would mimic the actual coupling found for the normally smaller densities.Indeed, the mean free path between ions and neutrals in the two-fluid reconnection simulations of Leake et al. (2013), which was self-consistently calculated, was about 100 m.This is nicely situated between our 2fl and 2flα case.Moreover, there are different recipes for calculating collisional frequencies, due to different values for the collisional cross sections available in the literature, which lead to large differences (Martínez-Sykora et al. 2017).The very small mean free path in our 2fl case suggests that it can be considered an MHD limit case, and this will be confirmed in our simulations below.

Results
In all six cases considered (MHD, 2fl, 2flα for the coronal versus TR reconnection case) the reconnection develops, producing bidirectional outflow jets, traveling away from the reconnection point, which can be seen in the snapshots of the total density shown in Figure 3.As the atmosphere is gravitationally stratified, the jets traveling upwards are denser and those traveling downwards are less dense than the surrounding fluid located at the same height.The figure shows that for TR reconnection, we get a pronounced upwards jet that is accompanied by a vertically oriented CS that lengthens as time progresses.From the time evolution of the upwards moving jets in the z rec = 2 Mm case we can estimate a vertical velocity of ≈ 8 km/s.An online animation for this case (z2.0-2flaplha.mp4)overplots the adaptive grid for the neutral density.In the case of coronal reconnection, we find a clear two-sided (up-down) jet forming, where the lower one ultimately interacts with the TR and chromosphere (forming post-flare loops).For both reconnection heights considered, z rec = 2 Mm and z rec = 4.5 Mm, the MHD (top row) and the 2fl models (middle row) give very similar results, as expected according to the collisionality regime.A clear difference appears in the snapshots for the 2flα model (bottom row of Figure 3), with the development of a less dense region with an edge of increased density, and this is seen in both z rec cases.In order to better understand these differences, we further analyze snapshots for z rec = 2 Mm at t = 622.6 s, the last time shown in Figure 3 (left column).

Analysis for the case with TR reconnection
We plot in Figure 4 the out-of-plane current density J y for the 2fl and 2flα models, and overplot magnetic field lines and isocontours of total density.Except for the fact that the 2flα seems a bit further evolved in time (as the top magnetic island is located a bit higher), the current density structures are very similar for the two models.We observe that the dense edges in the 2flα case are located just outside the current sheets.Because of the localized resistivity, the simulations evolve towards a Petschek-like reconnection.Slow shocks that accompany the reconnection, traveling in the x-direction, widen gradually the distance over In order to understand this clear difference in the 2flα case, which must be caused by two-fluid effects, we plot in Figure 7 different quantities along a horizontal cut located at z = 4 Mm.These 1D profiles are consistent with the previous 2D images, most notably Figures 5 and 6, where they cut across the V-shaped current sheet structures discussed earlier.For the 2fl model the structures in the charged and neutral density (panels (a) and (b)) are similar.The x-velocity (panel (c)) and temperature (panel (d)) profiles overlap for both fluids, meaning that the two fluids are coupled in both velocity and temperature.For the 2flα model, more charges imply less neutrals and higher density implies lower temperature.The temperature of the neutrals increases by more than 1.6 × 10 6 K at the center of the CS.Further from the center of the CS, the velocities of charges and neutrals are coupled, however inside the CS, they are completely different.The neutral velocity changes sign at the center of the CS, meaning that the neutrals are going out of the CS.Hence, the 2flα case which shows a pronounced anticorrelated structure in density and Article number, page 5 of 16 A&A proofs: manuscript no.main temperature for charges versus neutrals is clearly demonstrating decoupling between both species through the CS structure.
Panels (e) and (f) in Figure 7 show the total density, ρ, and the temperature of the center of mass, defined as The total density profile is very similar to the separate neutral and charged density profile for the 2fl model.For this case, this positive peak in the density, and a corresponding negative peak in the temperature, are related to the fact that the (vertically flowing) reconnection outflow comes from a lower height with higher density and lower temperature.However, the two-fluid solution in the 2flα case, Article number, page 7 of 16  where the collisional effects are enhanced, behaves rather different, and demonstrates a nonlinear runaway effect.In the case of the 2flα model, there is a central depletion in total density, where charges accumulate and neutrals deplete, and this is surrounded by a layer of enhanced total density where charges deplete and neutrals accumulate.A tiny central positive peak, of similar magnitude as in the 2fl model, still remains in the very center CS region.Seen from an MHD point of view, the whole fluid (both neutrals and charges) heats because of the collisions.This creates a peak in the center of mass temperature (panel (f), dotted line) and the entire fluid expands, so that the total density decreases towards the center of the CS (panel (e), dotted line).In the various panels of Fig. 7, vertical dotted lines show extremal positions in the neutral collisional heating term for case 2flα, discussed in more detail in what follows.
In order to better understand the temperature and velocity profiles seen in Fig. 7 we show the terms which enter the temperature (top panel) and velocity (bottom panel) equations in Fig. 8.The primary and secondary maximum peaks in the neutral frictional heating term (blue dashed line) for the 2flα case (top-right panel), are located at a distance of 0.03 Mm and 0.1 Mm away from the center of the CS, and these are indicated by vertical violet and gray dotted lines, respectively (those are repeated in each frame of Fig. 7).The peak in the Ohmic heating term (red solid line) is located close to the primary maximum peak in the neutral frictional heating term of the 2flα case (i.e.near the vertical violet dotted line).This frictional heating is negligible in the 2fl case (top-left panel).In the 2flα case (top-right panel), the charged fluid frictional heating term (red dashed line) has a similar profile to the neutral frictional heating term, but with smaller values, because of the higher charged density compared to the neutral density.In the 2flα case, the neutral frictional heating peak is five times larger than the peak value in the Ohmic heating.
We now turn attention to what causes the velocity decoupling effects, by discussing the forces as shown in Fig. 8, bottom panels.The inflow into the CS (as seen in panel (c) of Fig. 7) is driven by the magnetic pressure gradient which acts only on charges, producing initial decoupling in the velocities between neutrals and charges when entering the CS.This magnetic pressure gradient (red dashed line in bottom panels of Fig. 8) pushes the charges and the collisionally coupled neutrals inside the CS against the gradient of their pressures (solid lines).In the 2fl case, only Ohmic heating matters, and the overall decoupling between neutrals and charges stays minimal throughout the CS.In the 2flα case, the neutral and charged fluids decouple in velocity in a more pronounced way throughout the entire CS, again indirectly driven by the magnetic forces which drive the reconnection.But the initial decoupling grows and collisional heating of neutrals causes the neutrals to heat, expand and go out of the CS, accumulating between the primary and secondary heating point (violet and gray lines, respectively), as observed in panel (b) of Fig. 7.The inflow of charges there decreases because of the increased collisions with the neutrals that accumulated outside the CS (see the local minimum observed for the dotted red curve at the location of the gray line in panel (c) of Fig. 7).Therefore, the decoupling in velocity increases further, and also the frictional heating at this secondary neutral heating maximum point (the gray vertical line), which again increases the neutral pressure.Thus, the neutral inflow increases as a local maximum is observed for the dotted blue curve at the location of the gray line in panel (c) of Fig. 7.This is because of the neutral pressure gradient, which shows a local maximum peak in the solid blue curve at the location of the gray line in bottom-right panel of Fig. 8.This in turn drags the charges into the CS: a local maximum peak is seen for the red dotted line in bottom-right panel of Fig. 8, which has positive values around the gray line.This implies acceleration of the plasma towards the center of the CS, in the same direction and partially overlapping the magnetic pressure gradient curve.Therefore, the velocity of charges towards the center of the CS is increased and has  opposite sign to that of the neutrals, leading to more decoupling and to a runaway nonlinear instability.The heating of the charges at the primary heating point (violet line) slows down in this runaway instability process.The charges enter faster into the CS because of decreasing collisions with neutrals and this accelerates the process (besides leading to a thinner CS and faster outflows).The runaway process is further enhanced by the charges expanding due to the (collisional) heating at the secondary point (gray line).The runaway process creates a region of accumulation of neutrals bordering the CS (within the hydrodynamic timescale on which the charges enter in the CS due to the magnetic pressure gradient), which creates a secondary (collisional) heating point whereby the charges get pushed faster into the CS.

Analysis for coronal reconnection
We next study the case z rec = 4.5 Mm where the reconnection was triggered in the coronal region.A zoomed view shown in Figure 9 shows what happens near and above the TR, after the downwards reconnection outflow hits the dense material in the chromosphere and is reflected.The snapshots of density of neutrals versus charges (top row) show again a pattern of evacuated versus enhanced regions with dense versus evacuated edges.The decoupling in velocity (bottom left panel of Figure 9) points across the magnetic field lines and the isocontours of the decoupling in temperature (T c − T n ), shown in the bottom right panel, follow the structures in current density, suggesting that the decoupling is related to the magnetic fields.Smaller cur-rent sheets are formed at different locations and with different orientation, meaning that this separation of neutral versus charged densities across the current sheet is generally related to reconnection and not determined by gravity or localized resistivity.In a time evolution (an animation z4.5-2flalpha.mp4 is provided) we see that this process of separation seems to reverse in this case, when the structure rises again, but this is influenced by mixing with both neutral and charged material coming from below.Hence, in both 2flα cases, we find clear evidence of a nonlinear runaway process that enhances the spatial separation between charged and neutral fluids near current sheets.This runaway process is clearly related to having a large collisional free path, as compared to the width of the CS, as it is not seen at all in the 2fl model.Overall, we identified a runaway (decoupling) instability in a weakly collisional regime, which occurs in a nonstationary two-fluid setup to the reconnection problem.This will ultimately break down the fluid assumptions if neutral and charged densities decrease towards zero in separate locations.The exact conditions for the onset of this nonlinear instability, such as the α collisional parameter range, can be subject of a more idealized (without gravity) setup in the future.

Reconnection rate
We then calculate the reconnection rate Article number, page 9 of 16 A&A proofs: manuscript no.main in the same way as Leake et al. (2013); Murtas et al. (2021).The superscripts indicate the points where quantities are evaluated.The reconnection point X is defined as the point where the dissipation ηJ 2 y is maximum (we included η in this calculation because of the localized resistivity).X will be located at the center of the current sheet (x = 0) at some height, which might not be the initial reconnection point z rec , because this reconnection point migrates vertically, especially in the TR reconnection case z rec = 2 Mm, where the stratification is stronger.Initially, the X-point is located at the height z = z rec , which is one of the parameters of the simulations.The border of the CS, indicated by the superscript up , is defined as the point located at the same height as X, but displaced along the x-direction, as it is located at the point where the current density is half the current density measured at X.Because of the horizontal symmetry, the values calculated at the two opposing symmetric points at the borders of the CS are averaged, because numerically exact left-right symmetry might not be preserved.The Alfvén velocity v * A is calculated using the total density evaluated at the reconnection point X and the magnetic field B up z .We plot the thus computed reconnection rate as a function of time in Figure 10, for both TR and coronal reconnection, each time comparing MHD and both two-fluid rates.The reconnection rate has an initial increase, which is steeper for z rec = 4.5 Mm than z rec = 2.0 Mm.This is because at z rec = 4.5 Mm the density is lower and the CS thins faster, however the minimum width of the CS achieved during the simulations is similar for both heights.The maximum reconnection rate reached after this increase phase is M ≈ 0.1, a typical value for reconnection scenarios which use localized resistivity.Then, it decreases for both cases, mainly because the magnetic field dissipates.The reconnection rate decreases much faster in the TR reconnection case that started at z rec = 2.0 Mm, because of the quick and continuous upwards migration of the X-point over a distance of ≈0.5 Mm at the end of the simulation.It is seen that the reconnection rates remain similar for the MHD and the two two-fluid setups.It is known that the plasmoid formation increases the reconnection rate (Leake et al. 2012;Murtas et al. 2021), however, in our case the plasmoids appear too late in the simulations to influence significantly the growth rate.

Secondary plasmoids
In the last snapshots of the simulations with TR reconnection (z rec = 2 Mm) we can observe the formation of secondary plasmoids, but not for the coronal case z rec = 4.5 Mm.Their shape can be clearly seen in the current density map in Figure 4 for the 2flα case, where we can visually estimate a size of ≈0.3 Mm.The initial formation phase of the plasmoids can be best seen in the vertical profile of the vertical velocity.Several equidistant moments of time are shown in Figure 11 for the three models: 2flα, 2fl and MHD.In these profiles, the plasmoids are seen as secondary velocity extrema that develop and get advected upwards.The plasmoids move with the outflow from the primary reconnection point, which are at the locations where the vertical velocity changes sign (fitted locally by linear slope in the three panels).We can visually estimate that the fastest growing length scales are similar in the three cases, however the growth is largest in the MHD case, followed by 2flα and then 2fl.This suggests that the two-fluid effects do not affect the initial phase of the growth of these plasmoids.We do expect that the later stages of this plasmoid evolution will be similarly affected by the runaway decoupling process in the 2flα case, when charges accumulate towards the middle of the CS.
The linear growth of the tearing mode is affected by the gradient in the vertical flow (parallel to the CS), as shown in a linear resistive MHD analysis about a 2D configuration by Bulanov et al. (1978).In their analysis the vertical profile is assumed linear v z (z) = az, where a quantifies the vertical velocity variation.Bulanov et al. (1978) show that the growth rate γ(a) of the tearing mode at a finite a, is reduced from the case without flow gradient, a = 0, in the sense that γ(a) ≈ γ(0) − a.
A linear analysis of the tearing mode in the two-fluid approach, using simplified assumptions of uniform density, and hence no gravity (and no added vertical flow) is presented in Appendix A. This shows that two-fluid effects will define an effective density ρ c0 ≤ ρ eff ≤ ρ T0 between the background charged ρ c0 and total ρ T0 densities.The linear growth of the tearing mode in this simplified twofluid assumption is bounded between the growth calculated in the MHD assumption when the density is the total density (ρ T0 ) and the charged density (ρ c0 ).In our setup the Article number, page 10 of 16 B. Popescu Braileanu and R. Keppens : Two-fluid reconnection jets in a gravitationally stratified atmosphere.ionization fraction is high, and the difference that the twofluid effects might introduce in the growth of the tearing mode are of the order of 10 −4 s −1 .This is shown in our appendix, and can be seen as the difference on the y-axis of Fig. A.1 between the orange point which corresponds to the collisional parameter used in 2flα simulations and the red point which correspond to a value of the collisional parameter larger than the maximum value of α in 2fl cases.This two-fluid related difference is one order of magnitude smaller than the obtained difference in the velocity slope, which we show for the three cases in Figure 11.The ordering of the growth rates of the plasmoids, as visually estimated from Figure 11, is reversed compared to the ordering of these slopes.This is consistent with the expected reduction in tearing growth rates due to added vertical velocity gradients.Similarly, the large (vertical) gradients in the outflow Article number, page 11 of 16    velocities are probably the reason why we do not observe secondary plasmoids in the simulations with z rec = 4.5 Mm, since then the vertical velocities and corresponding vertical gradients are larger than for the z rec = 2 Mm case.Therefore, in our simulations, the initial linear growth of the plasmoids is rather influenced by the fact that they form at different moments and slightly different heights during the simulations, when the vertical gradients in the velocity are different.The difference in the growth rate produced by this gradient in velocity is much larger than the difference in the growth introduced by the collisional effects.

Synthetic views
As a relatively straightforward observational validation, we can produce synthetic images resulting from emission from optically thin spectral lines.To do so, we calculate the emission in the Solar Dynamics Observatory Atmospheric Imaging Assembly (SDO/AIA) channel AIA 193 Å, which has its highest response for temperatures between 10 6 and 2 × 10 6 K.We do so only above the height z ≈ 2 Mm in our initially stratified atmosphere, because the chromospheric region itself is not appropriately dealt with in an optically thin limit.Being optically thin emission, we deal essentially with a function of local temperature T and number density n given by where R(T ) is a response function which depends on temperature, and a log R − log T view is shown in Figure 12 for the wavelength 193 Å, obtained using the CHIANTI atomic database (Del Zanna et al. 2015).Images in these Extreme Ultra-Violet (EUV) channels of AIA can be synthesized from 3D data cubes by integrating Eq. 12 along the line of sight (LOS).We have 2.5D data, with an assumed invariance in the third (y) direction.As we have 2.5 D simulations, the images are shown at the simulation resolution, and there is no integration along LOS, also making units meaningless.As the AIA 193 channel is emitted by Fe XII, emission is obviously from charged species.From our twofluid plasma-neutral model, we must choose which number density and temperature is taken in Eq. 12, and a more accurate result might be obtained by using only the charged density, instead of the total density.In particular, plasmoids are usually observed as bright blobs due to their increased density compared to the surrounding medium.In this case models which consider the total density in the calculation of the synthetic image (as is usually done from a single-fluid MHD simulation) would give wrong results, opposed to considering the charged particle density.Figure 13 shows synthetic images in the AIA 193 channel which capture the temporal evolution of the plasmoids.We found secondary plasmoids when z rec = 2 Mm, and show synthetic views for the 2fl case (top row) and the 2flα case (middle row), where we use the charged fluid temperature and density, meaning adopt Λ(n c , T c ) in Eq. 12.The bottom row shows for the case 2flα the same snapshots as in the middle row, but instead using the total density and the center of mass temperature defined by Eq. 10, in practice using Λ(n T , T 2fl ) in Eq. 12.
On all panels of Figure 13, we show an emission isocontour at a fixed value.We can observe that the plasmoid fades away as the surface delimited by the isocontour becomes smaller, while it accelerates upwards in the 2fl case (top row).On the contrary, the plasmoid decelerates and becomes brighter in the 2flα case (middle row), most probably because of the increasing charged density due to the runaway effect.A completely different interpretation results from the images constructed using total density for the 2flα case (bottom row), where the plasmoids appear as dark structures that get surrounded by two bright and descending spikes, coming down from z ≈ 4 Mm, where the accumulation of the neutrals outside the CS seems to be maximal (Figure 4).Because of the runaway process, the synthetic image which uses total density gives a wrong image, as the neutrals accumulated outside the current sheet (considered in the total density) would not actually contribute to the emission in this line.
Therefore, charged density only should be used in generating synthetic images in optically thin lines in order to produce more accurate results.Finally, we investigate how the MHD model estimates the charged density, and quantify differences obtained in synthetic images for the three models.The ionization fraction depends on height only and is constant in the MHD model, depending on the background profile.The MHD model cannot track how the ionization fraction changes during the simulation because of how the two species evolve, but we can define a quantity R, being the inverse of the normalized mean molecular weight, as calculated using the two-fluid equilibrium (hence the 0 subscripts) atmosphere which can be used to retrieve the temperature in the MHD model in a consistent way with the two-fluid model As we also know the mean molecular weight of charges and neutrals at each height1 , then we can estimate the charged and neutral density from the MHD model: With these identifications, we can mimic consistent twofluid like quantities from a pure MHD model, and then also turn the MHD in a synthetic image based on temperature and charged density.In practice, for the MHD model the density of charges is then estimated using Eq. 15. Figure 14 shows the resulting images for the three models and for the two cases of reconnection heights.We can observe that the MHD and 2fl models give very similar results.There is a small difference: in the z rec = 2 Mm the 2fl shows less intensity than the MHD, but the reverse happens for the z rec = 4.5 Mm case.This is due to the fact that the density of charges is slightly overestimated for the fluid coming from below and underestimated for the fluid coming from above.Because of different scale heights between charges and neutrals, the ionization fraction is different at different heights.This is an intrinsic limitation of the MHD model which only keeps the information of the initial ionization state.The snapshots used for z rec = 4.5 Mm are at an earlier time than the final time of this simulation, where we see a clear upward jet feature as the reconnection outflows impact the chromosphere.This jet-like feature with width comparable to the width of the CS (≈ 30 km) is present in all three images, and may be observable in high cadence, high resolution observations.

Summary
We did simulations of two-fluid reconnection in a gravitationally stratified atmosphere, similar to the solar atmosphere, where we studied the collisional effects for reconnection points situated at different heights.Because of the localized resitivity used, a Petschek type reconnection developed with slow shocks propagating along the x-direction, disrupting the current sheet in a V shape.The vertical upward velocities of ≈ 10 km/s in the z rec = 2 Mm and the width of the jets, which vary from being similar to the width of the CS (≈ 30 km) near the X-point to ≈ 700 km, higher up, due to the slow shocks, are consistent with properties of type I spicules.The MHD and 2fl simulations showed very similar results.When z rec = 4.5 Mm the reconnection outflow hit the denser material below and interacted with reconnected magnetic field, creating secondary thin current sheets, leading to locally more turbulent behavior in the post-flare loop region.
The thermal effect of the collisions on the evolution of the neutral fluid has been observed earlier in simplified 1D slow shock simulations, where the heating of the neutrals produces an overshoot in the neutral velocity (Hillier et al. 2016;Snow & Hillier 2019).In our case, the decoupling is larger when the collisional effects are increased (the 2flα cases) and the heating of the neutrals will produce a runaway effect which separates the neutrals and the charges across the CS.The neutrals accumulate outside the CS, while the charges tend towards the center of the CS.This gives reversed contrasts in the charged density and neutral or total density maps.The regions with very low density of neutrals have high density of charges (inside the CS) and the regions with very high density of charges have low density of neutrals (outside the CS), and these differences increase over time.The temperature maps are consistent with the density maps, with high density regions having low temperatures, while low density has high temperatures.This was analyzed in detail, and a nonlinear decoupling runaway effect was identified.
We obtain high reconnection rates in all the simulations because of the localized resistivity.The localized resistivity has the effect to bound the current density, similar to an anomalous prescription.We obtain the maximum reconnection rate in the Petschek model, which is 0.1 (Vasyliunas 1975).Two-fluid effects do not increase our reconnection rates further, as opposed to the results of Leake et al. (2012Leake et al. ( , 2013)).Our setup is closer to conditions relevant for the stratified solar atmosphere.
At later time we observe the formation of secondary plasmoids.They are observed for all the models when z rec = 2 Mm.The large outflows and associated gradients in the flow when z rec = 4.5 Mm are likely inhibiting the linear tearing mode.This effect is much more important than collisions in the early formation of plasmoids.In simplified assumptions of uniform density, collisions define an effective density in the two-fluid model, similar to the idea presented by Murtas et al. (2021) in simulations of the coalescence instability.At later stages, however, the accumula-tion of plasma adds to the nonlinear evolution of the tearing mode.Because of a smaller effective density when the collisional coupling is reduced, the effective Alfvén speed is larger and the simulations evolve slightly faster.
The secondary plasmoid formation process would look completely different in synthetic images which use total density instead of charged density in the calculation of emission in optically thin lines.When the charged density is used for the synthetic image generation, an MHD model which does not keep track of the ionization fraction produces slightly different results because of this limitation.It became evident that using charged particle densities leads to more realistic behavior, in line with observed enhanced emission blobs.Article number, page 16 of 16

Fig. 1 .
Fig. 1.Initial profiles as a function of height.Densities: neutrals (blue dashed line), charges (red dotted line) and total (black solid line) on the left axis.Temperature (solid orange line) on the right axis.
Fig.2.Collisional mean free paths as functions of height.We show the mean free paths between ions and neutrals (λin, red lines) and between neutrals and ions (λni, blue lines) for the two collisionality regimes considered: self-consistent from plasmaneutral parameters (2fl, solid lines, left axis) and reduced case 2flα (dashed and dotted lines, right axis, note the different order of magnitude).The collisional mean free paths are calculated according to Eq. (8).
Braileanu and R. Keppens : Two-fluid reconnection jets in a gravitationally stratified atmosphere.

Fig. 7 .
Fig. 7. Horizontal cut at coronal height z = 4 Mm at time t = 622.6 s for the simulations with TR reconnection or zrec = 2 Mm.(a) Density of charges, (b) Density of neutrals, (c) x-velocity for charges and neutrals, (d) Temperature of charges and neutrals, (e) Total density, (f) Center of mass temperature from Eq. (10).Different curves are for the two models considered: 2fl and 2flα and quantities for charges and neutrals are shown in the same plot in panels (c) and (d), as indicated in the legend.The primary and secondary maximum peaks in the neutral collisional heating term for case 2flα, as shown in detail in Fig. 8, located at z = −0.03Mm and z = −0.1 Mm are indicated by vertical violet and gray dotted lines, respectively.

Fig. 8 .
Fig. 8. Horizontal cuts along z = 4 Mm for the case with TR reconnection or zrec = 2 Mm at time t = 622.6 s, the same time as for Figure 7. Quantities corresponding to charged fluid are shown with red lines and those for neutral fluid with blue lines.Left panels: 2fl case; Right panels: 2flα case.Top row: Terms in temperature equation: Ohmic term (solid red line), collisional heating terms (blue dashed lines: red for charges and blue for neutrals).Bottom row: Terms in velocity equation; acceleration due to: magnetic pressure (red dashed line), gas pressure gradient (solid lines: red for charges and blue for neutrals), collisions (dotted lines: red for charges and blue for neutrals).The primary and secondary maximum peaks in the collisional heating associated with neutrals, located (and shown only for negative x, because of symmetry) at z = −0.03Mm and z = −0.1 Mm are indicated by vertical violet and gray dotted lines, respectively for the 2flα case (right panels).

Fig. 9 .
Fig. 9. Two-fluid coronal reconnection downflows impacting the chromosphere.Snapshots of charged density (top-left), neutral density (top-right), out of plane current density (bottom-left) and charged temperature (bottom right) for the simulation zrec = 4.5 Mm, 2flα at t = 591.5 s, the last time shown in the right column, bottom row of Figure 3.In the panels of the neutral and charged densities (top panels) the orange arrows show the velocity of neutral and charged fluid, respectively.We overplot the magnetic field lines and the decoupling velocities over the current density map (bottom-left panel) and isocontours of the temperature decoupling (Tc − Tn) in the range −1.8 × 10 6 K and 1.8 × 10 6 K over the charged temperature image (bottom right panel).

Fig. 10 .
Fig. 10.Reconnection rates for the six simulations.For zrec = 2 Mm (left panel) and zrec = 4.5 Mm (right panel) as a function of time.Different curves indicate different models: MHD (orange solid line), 2fl (green dashed line) and 2flα (blue dotted line).

Fig. 11 .
Fig. 11.Vertical profiles of the vertical velocity when the plasmoids form for zrec = 2 Mm.Shown for five moments of time in the late stage of the simulations, for the three models: 2flα (left panel), 2fl (middle panel) and MHD (right panel).For the two-fluid models (left and middle panels) we show the velocity of charges.The linear fit of the velocity profile around the reconnection point (where the velocity changes sign) is shown by a black line and the value of the slope is indicated at the top of each panel.

Fig. 13 .Fig. 14 .
Fig.13.Synthetic views on secondary plasmoids, at three consecutive times.For plasmoids in the case zrec = 2 Mm, we show synthetic images of AIA 193Å emission.Top row: 2fl case; middle and bottom rows: 2flα case.For the calculation of emission, we used the density and temperature of charges in the top and middle row, and the total density and the center of mass temperature of neutrals and charges for the bottom row.The units are non-dimensional.We kept the colorbar in the figures, for comparison purposes, as normalization is done in the same way.The last row shows different limits on the colorbar because of the increased density when total density is considered instead of the charged particles only.