Fast magnetic reconnection plays a fundamental role in driving explosive dynamics and heating in the solar chromosphere. The reconnection time scale of traditional models is shortened at the onset of the coalescence instability, which forms a turbulent reconnecting current sheet through plasmoid interaction. In this work, we aim to investigate the role of partial ionization in the development of fast reconnection through the study of the coalescence instability of plasmoids. Unlike the processes occurring in fully ionized coronal plasmas, relatively little is known about how fast reconnection develops in partially ionized plasmas (PIPs) of the chromosphere. We present 2.5D numerical simulations of coalescing plasmoids in a single fluid magnetohydrodynamic (MHD) model and a two-fluid model of a partially ionized plasma (PIP). We find that in the PIP model, which has the same total density as the MHD model but an initial plasma density two orders of magnitude smaller, plasmoid coalescence is faster than the MHD case, following the faster thinning of the current sheet and secondary plasmoid dynamics. Secondary plasmoids form in the PIP model where the effective Lundquist number $S=7.8\xd7103$, but are absent from the MHD case where $S=9.7\xd7103$: these are responsible for a more violent reconnection. Secondary plasmoids also form in linearly stable conditions as a consequence of the nonlinear dynamics of the neutrals in the inflow. In the light of these results, we can affirm that two-fluid effects play a major role in the processes occurring in the solar chromosphere.

## I. INTRODUCTION

Magnetic reconnection is a process responsible for explosive events in astrophysical, space, and laboratory plasmas, allowing plasmas to break free from the frozen-in constraint imposed by ideal magnetohydrodynamics (MHD).^{1} The magnetic field lines reconnect through a narrow diffusion region, enabling the conversion of stored magnetic energy into kinetic and thermal energy, and particle acceleration.^{1,2}

A classical description for reconnection is provided by the Sweet–Parker model,^{3,4} which describes steady-state reconnection in a long diffusion region of length *L* and width *δ*. The reconnection time scale goes as $\tau SP\u223cL/\delta \u223cS1/2$, where $S=LvA/\eta $ is the Lundquist number, *η* is the diffusivity, and *v _{A}* is the Alfvén speed.

The presence of localized, transient outflows in the solar chromosphere (such as the chromospheric jets^{5,6}) is believed to be the result of the onset of magnetic reconnection. Further evidence for magnetic reconnection is provided by the identification of bubbles of plasma in the outflow of chromospheric jets and UV bursts,^{7–9} which are generally interpreted to be plasmoids.

Plasmoids are commonly present in reconnecting systems,^{10–16} and they are believed to play a major role in speeding up reconnection by having an influence on the variation of the current sheet size.^{17} Under the formation of plasmoids, the current sheet breaks into fragments, and the resulting high current densities in each of these sections facilitate a high reconnection rate.^{18} Plasmoid formation due to the instability of Sweet–Parker current sheets has been extensively examined through numerical studies^{16,19–25} and found in direct imaging observations of solar flares.^{26} Several works proved that in fully ionized plasmas it is possible to trigger plasmoid formation in thin current sheets ($\delta /L\u226a1$) for a critical Lundquist number $\u223c(103\u2212104)$.^{13,15,16,18,21,25,27–32} Numerical studies^{27} found an upper limit for the critical Lundquist number $\u223c3\xd7104$. Above this critical value of *S*, the current sheet becomes rapidly unstable to the resistive tearing instability, forming a chain of plasmoids along the current sheet.^{14}

In astrophysical plasmas with very large *S*, the onset of the tearing mode takes place in current sheets that collapse to a thickness on the order of $S\u22121/3$, larger than the Sweet–Parker thickness of $S\u22121/2$. The trigger of fast reconnection before a Sweet–Parker-type configuration can form was examined in detail by several studies.^{33–36}

Many works proved that the critical Lundquist number and aspect ratio for the onset of the tearing instability are affected by changes in the initial setup, which might result in a discrepancy of orders of magnitude in the critical Lundquist number.^{32,37,38} A major role is played by the initial current sheet configuration^{31,32,37} as well as by the amplitudes of viscosity and perturbation noises,^{37,39–41} and the plasma *β.*^{31,32,37,38} The variation of the critical *S* as a function of the initial noise, investigated in some works,^{37} covers several orders of magnitude (from $S\u223c103$ to $S\u223c106$). The role of the configuration of the simulation domain in affecting the threshold of the critical Lundquist number has been pointed out in studies of magnetic reconnection in fully ionized plasmas.^{31,32} In these works, it was shown that a 2.5D simulation of magnetic reconnection with a force-free current sheet and uniform plasma pressure as initial conditions lead to a much lower critical *S* than those obtained in 2D cases with an initial Harris current sheet and nonuniform plasma pressure. Their results are consistent with the findings of a recent study,^{38} where magnetic reconnection is examined on different spatial scales in weakly ionized plasmas by using a reactive 2.5D multi-fluid plasma–neutral model.^{42,43}

The nonlinear tearing mode shows the development of an important secondary instability called the coalescence instability.^{2} This instability is driven by the coalescence of neighboring plasmoids sharing an X-point and results from the attractive forces between parallel currents. The coalescence instability is characterized by two different phases: ideal and resistive. The ideal phase has a growth rate that is almost independent of *η.*^{44} The resistive phase is driven by the current sheet reconnection.^{2} In a single-fluid MHD approach, the coalescence instability produces a reconnecting current sheet by driving plasmoid interaction and is a key process that might explain fast reconnection without the need of anomalous resistivity terms to be added into the system.^{18}

In the solar chromosphere plasmas, as well as many other plasmas found in the universe, are partially ionized, and their ionization degree falls in the range of $10\u22124\u221210\u22121$.^{18,45–47} Multi-fluid effects linked to the different behavior of the particle species must be taken into account for a correct physical description of this atmospheric layer. The low chromospheric densities do not allow a complete collisional coupling between ions and neutral species: the low ion fraction allows the fewer charged particles to be coupled to the neutrals, but the neutrals may not be completely coupled to the ions. A partial coupling between the two species results in the presence of relative motions. While ideal MHD equations are applicable to fully ionized plasmas, two-fluid effects should be described by taking into account the relative motions between plasma and neutral components in the solar chromosphere. The role of partial ionization in the onset of magnetic reconnection and development of the resistive tearing instability was investigated in many studies.^{7,17,18,42,43,48,49} In a system where the two fluids are coupled through elastic collisions and subject to ionization and recombination, the reconnection rate in the coalescence process depends on the ion fraction.^{50,51}

The role of fast magnetic reconnection in triggering events in the solar chromosphere is undoubtedly fundamental. An additional complexity comes from the partially ionized nature of the solar chromosphere. In this paper, we discuss the role of partial ionization in magnetic reconnection through the study of plasmoid coalescence, with the aim of understanding to what extent the two fluid effects influence such process. In order to do this, we first compare two reference MHD and partially ionized plasma (PIP) simulations (Sec. III) and then investigate in more detail how two-fluid properties affect the coalescence instability through a parameter survey (Sec. IV). In Sec. V, our results are connected back to the physical scales of reconnection in the solar chromosphere.

## II. METHODS

We perform simulations of the coalescence instability in fully (MHD) and partially ionized plasmas (PIP), using the (PIP) code.^{52} The code makes use of a four-step Runge-Kutta and a fourth-order central difference scheme. The fully ionized plasma consists of a single-fluid model of a hydrogen plasma. The partially ionized plasma environment is simulated through a two-fluid model consisting of two separate sets of equations, describing a neutral fluid and a charge-neutral ion–electron plasma which are collisionally coupled. The equations are derived from those found in previous models.^{42,53,54}

All sets of equations are non-dimensionalized. The choice of this particular normalization is performed in order to have a dependency on characteristic length scales, which are comparable to the size of the plasmoids involved in the merging. This allows the model to be applied to plasmoids at different scales in the solar chromosphere, from a few meters to a few hundred kilometers, as all the quantities scale together with the plasmoid size. Such normalization also affects the collisional coupling between the two fluids. As a characteristic dimensional time scale *τ* can be derived from the physical plasmoid size and the characteristic speed in the solar chromosphere (in our case this is the sound speed), the non-dimensional collisional frequency is easily re-scaled back to physical quantities by dividing it by *τ*. Further details on the non-dimensionalization are provided at the end of this section and in Sec. II A.

The neutral fluid is described by non-dimensional inviscid hydrodynamics equations

while inviscid resistive magnetohydrodynamics relations govern the plasma, here stated in non-dimensional form

In the equations above, the subscripts *p* and *n* refer, respectively, to the ion–electron plasma and the neutral fluid, **v**, *p*, *ρ*, *T*, and *e* are the velocity, gas pressure, density, temperature, and internal energy of each species, $\gamma =5/3$ is the adiabatic index, and **B** is the magnetic field. Both fluids follow the ideal gas law. The factor 2 in Eq. (12) is to take into account the electron pressure. The parameter *α _{c}*, given by Eq. (13), is associated with the two fluids collisional coupling. This particular formulation of

*α*is new to models of magnetic reconnection in partially ionized plasmas, and it accounts for the increased amount of collisions at supersonic drift velocities. The non-dimensional expression for

_{c}*α*

_{c}^{55}including charge exchange

^{56}is as follows:

where $\alpha c(0)$ is the initial coupling and *v _{D}* = |

**v**-

_{n}**v**| is the magnitude of the drift velocity between the neutral components and the hydrogen plasma. When the drift velocity becomes bigger than the thermal velocity, the particles are subject to a higher number of collisions as they are drifting past each other. The collisional coupling between ions and electrons is represented by setting a small finite diffusivity

_{p}*η*that is assumed to be spatially uniform and not varying with time. In this work, we are not including the Hall effect.

The two systems of equations are non-dimensionalized^{57} by a reference density *ρ*_{0} and the total sound speed $cs=\gamma (pn+pp)/(\rho n+\rho p)$, initially set equal to 1. For the MHD simulation, where the plasma is fully ionized, the initial density and pressure are constant and equal to

For the PIP simulations, the bulk density and pressure are equal to the MHD values

and they are uniform in all the domain. Initially the two fluids are in thermal equilibrium.

### A. Initial conditions

The bulk initial conditions of the PIP case are equal to the initial conditions of the MHD case. The initial setup is provided by a force-free modified Fadeev equilibrium.^{2,58} The magnetic scalar potential of the classical 2D Fadeev equilibrium in the $xy\u2212$ plane is given by^{2}

where $B\u221e$ is the field intensity for the limit $|y|\u2192\u221e$. In our simulations, $B\u221e$ is equal to $2\gamma \u22121\beta \u22121$, where plasma $\beta =2p/B2,\u2009k=\pi 2$, and we set $\u03f5=0.5$, which corresponds to a moderately peaked current localization at the plasmoid center, shown in Fig. 1. As $\u03f5\u21920$ there is a weaker localization and a weaker attraction between the plasmoids, while $\u03f5\u21921$ corresponds to a peaked localization and stronger attraction forces. At the upper limit (*ϵ* = 1), the current distribution becomes the delta function.

In the mid to upper chromosphere the plasma *β* may become very small. Although the photospheric magnetic field emerging from the convection zone is not force-free, its structure is rearranged by the time it reaches the corona as the non-force-free components decay due to the action of chromospheric neutrals.^{59} It is hence of interest to study the coalescence instability in a regime that is initially force-free. The magnetic field *B _{x}* and

*B*components from the classic Fadeev equilibrium are not sufficient to satisfy the condition $J\xd7B=0$ for a force-free field. Therefore, we modify the traditional Fadeev equilibrium by including a component

_{y}*B*, making it force-free. The magnetic field components are shown in the following equations:

_{z}Setting *ϵ* = 0 in the equations above, the *B _{y}* component leads to a current sheet with the characteristic tanh profile of the well known Harris sheet.

^{60}Our initial conditions for the current density are displayed in Fig. 1 and are the same for both MHD and PIP cases.

The Fadeev equilibrium is unstable to the coalescence instability.^{2} We hence choose a small perturbation in the velocity of both plasma and neutral components to break the initial equilibrium by pushing neighboring plasmoids toward each other. The velocity perturbation is given by

where $vnoise$ is a white noise component simulating small environmental perturbations. The sine term dependent on *x* in the main perturbation results in a push on the pair of plasmoids, so they move closer to each other. As the domain has periodic boundary conditions on the sides, there is an effective chain of plasmoids moving along the *x*-axis. The perturbation causes the coalescence to take place for each pair of plasmoids separately, while moving the other plasmoids away. The term dependent on *y* localizes the perturbation to a small region around the plasmoids center. The white noise perturbation, which is set equal for plasma and neutrals, has a magnitude of 0.0005, two orders of magnitude smaller than the main perturbation in Eq. (22). Choosing such value prevents the noise from dominating the motion of the two plasmoids during coalescence, but allows the development of dynamics at a smaller scale by breaking the symmetry of the system. The same random noise seed was used in all simulations.

The reference MHD simulation in Sec. III is resolved by 2062 × 3086 grid cells, corresponding to a cell size of $\Delta x=1.95\xd710\u22123$ and $\Delta y=2.6\xd710\u22123$. In the PIP case, the partial ionization effects lead to the thinning of the current sheet and the development of sharp small-scale magnetic structures as a result of the neutrals decoupling from ions and leaving the current sheet.^{48,59,61,62} Our simulations show the formation of a thinner current sheet in between the plasmoids merging due to this two-fluid effect. The reference PIP case in Sec. III was hence run at the higher resolution of $\Delta x=1.2\xd710\u22123$ and $\Delta y=1.6\xd710\u22123$ to ensure the current sheet is resolved by our grid. The grid in the PIP case is composed of 6478 points in the *x* direction and 4862 points in the *y* direction. In Sec. IV, we present a parameter study of the coalescence process. The resolution used for each simulation and the total number of grid points are detailed in Sec. IV.

The initial separation between the plasmoids (calculated from $O\u2212$ point to $O\u2212$ point, identified as blue spots in Fig. 1) is equal to 4*L*, where *L* is resolved by 515 grid points in the MHD case and by 809 grid points in the PIP case. The plasmoid width, calculated as the distance between top and bottom edges of the separatrix and which initial value is $1.66L$ (resolved, respectively, by 638 grid points in the MHD case and by 1037 grid points in the PIP case), is determined by the conditions of the Fadeev equilibrium for the magnetic field.

The non-dimensional diffusion length scale is calculated as $Ldiff=4\eta \tau $ for both the reference cases and the simulations of the parameter survey. Taking *τ* = 1, the diffusion length scale for the MHD and PIP reference cases is $4.5\xd710\u22122$ for a diffusivity $\eta =5\xd710\u22124$, while for the cases in the parameter survey, which are characterized by $\eta =1.5\xd710\u22123,\u2009Ldiff=8\xd710\u22122$. The approximate number of grid cells per diffusion length scale at the lower resolution of the MHD case is (23, 17) in (*x*, *y*). These, respectively, increase to 41 grid cells along *x* and 31 grid cells along *y* for $Ldiff$ in the parameter survey. Therefore, the diffusion scale is resolved in all simulations. The collisional ion–neutral time scale $\tau col,pn=(\alpha c\rho n)\u22121$ for the reference PIP case is $10\u22122$, while the collisional neutral–ion time scale, defined as $\tau col,np=(\alpha c\rho p)\u22121$, is 1. These values lead to the non-dimensional coupling length scales $Lpn=10\u22122$ and *L _{np}* = 1, which are both well resolved by our grid.

### B. Boundary conditions

While coalescing, the plasmoids move toward each other along the *x*-axis. Because of the symmetry of the problem, in the reference MHD simulation (Sec. III) and in the set of simulations performed for the parameter survey (Sec. IV), we cut the computational domain at *x* = 0 and use symmetric boundaries, where *v _{x}* and

*B*change sign across each boundary and

_{y}*v*and

_{y}*B*remain the same. The left boundary is shown in Fig. 1 as a dashed green line. The computational domain size is chosen equal to $x=[0,4]$ and $y=[\u22124,4]$.

_{x}The dynamics of the plasmoids merging in the reference PIP case (Sec. III) is evaluated in a full computational domain, with $x=[\u22124,4]$ and $y=[\u22124,4]$. This arrangement was made to be able to better examine the dynamics in the region of the current sheet at higher resolution. In this case, the top and bottom boundaries are kept symmetric, while the side boundaries are chosen to be periodic.

## III. RESULTS

First, we explore the coalescence instability in both a single-fluid fully ionized plasma and a two-fluid partially ionized plasma by comparing two simulations, an MHD case and a PIP case. The single-fluid case acts as a reference case for the more complex two-fluid simulation.

The initial parameters for both simulations are constant diffusivity $\eta =0.0005$ and plasma $\beta =0.1$. In the PIP simulation, we set the collisional coefficient $\alpha c=100$ and the ion fraction $\xi p=\rho p/(\rho p+\rho n)=0.01$, while the effective ion fraction in the MHD case is $\xi p=1$. The effects of the parameters variation on the coalescence in PIP simulations are investigated later on in Sec. IV.

Figure 2 displays a sequence of the evolution of the current density *J _{z}*, which is directed out of the plane. For a better comparison, the frames show times where similar physics takes place in both the MHD and PIP cases. As $\eta \u22600$, the magnetic flux reconnects and leaves the current sheet that is formed in between the coalescing plasmoids: this is the region of strong negative current between the two plasmoids in panels (a), (b), (c) and (d) of Fig. 2. In a first phase, the current sheet length rapidly increases when the plasmoids approach, and then it progressively reduces with the size of the coalescing plasmoids [panels (c) and (d) of Fig. 2]. The reconnection results in the formation of a single large plasmoid, as shown in panels (e) and (f) of Fig. 2.

The left–right symmetry in the PIP case is broken during the coalescence, as evident, in particular, from panels (d) and (f) of Fig. 2. The asymmetry arises from the initial noise perturbation in Eq. (22), allowing the onset of small-scale dynamics in the central current sheet. The symmetry in the MHD case is also reinforced from the presence of a central boundary at *x* = 0, as introduced in Sec. II B.

Figure 3 displays in blue the separation of the two plasmoids in the MHD case, calculated as the distance between $O\u2212$ points. The squares along the curve identify the times of panels (a), (c), and (e) in Fig. 2. This distance fluctuates in time, with peaks that appear regularly during the reconnection phase. The merging plasmoids accelerate toward each other, move slightly apart as they bounce on the current sheet, and accelerate back again toward the center. Such movement can be associated with the high gas pressure generated inside the current sheet.

As observed from the black curve in Fig. 3, the distance between the two plasmoids reduces rapidly in the PIP case, following the faster reconnection process and without displaying the same oscillations that are remarkable in the MHD case. The reason could be that the reconnection is happening fast enough that the plasmoids do not need to rebound off each other. The squares here identify the times of panels (b), (d), and (f) in Fig. 2.

The shortening of the coalescence time scale in the PIP case can be associated with the decoupling of ions and neutrals in the reconnection region. This results in a faster thinning of the current sheet as the lower ion density allows a stronger compression than in the MHD case. The decoupling is discussed in more detail in Sec. III A. The sharpening of the magnetic field profile has already been observed in many studies^{61–66} as a result of the ambipolar diffusion. The reduced single-fluid model for two-fluid effects provided by the ambipolar diffusion provides a good approximation for strongly coupled systems, and it is consistent with the effects that we record in our simulations. However, it might fail in describing the complexity of weakly and intermediate coupled systems, and could not be used to explain the different time scale of the first phase of coalescence when the plasmoids move toward each other. The increased complexity of our system is, therefore, better investigated through a full two-fluid model such as the one used in this work.

At a first qualitative view, several differences are present between the fully ionized and the partially ionized cases. First, the plasmoid merging occurs faster in the PIP simulation, ending at $t\u223c4$, while coalescence takes a longer time in the MHD case, ending at $t\u223c18.5$. The end is identified with the time when a single, large plasmoid is fully formed and the current density at its center reaches a positive maximum, stabilizing to a constant value. The difference in the coalescence time scale is related to the reconnection rate, discussed in Sec. III B. A second difference is that there is no clear sign of shocks in the PIP case, while the MHD case shows an abundance of fine structures in panels (c) and (e) of Fig. 2. The identification and classification of these structures as shocks are investigated in detail in Sec. III C.

Two unique features are present in the PIP simulation only. The first is the production of secondary plasmoids in the central current sheet, linked to the onset of additional instabilities that will be discussed in Sec. III D. The second is the formation of an extended neutral jet in the reconnection region, whose properties are examined in Sec. III E.

### A. Ideal phase of coalescence

The mutual attraction of parallel currents pulls the two initial plasmoids together, and a current sheet forms as a result of the anti-parallel magnetic field being pushed together. In the MHD simulation, the ideal phase of coalescence is characterized by a plasma inflow forming along the *x*-axis that contributes to the formation of the central current sheet. The charged species in the PIP simulation are pulled to the center of the domain by the Lorentz force, and the neutrals are dragged by collisions resulting in a small drift velocity that can be seen in Fig. 4 for *t* = 1.62. The drift velocity ($vn\u2212vp$) differs from zero in the inflow, which indicates that the two species are weakly coupled in this first phase of coalescence, and increases steadily in time with the acceleration of the plasmoids motion toward each other. After *t* = 2.42, the reconnection process changes as the tearing instability takes place with the formation of secondary plasmoids [panels (b) and (d) of Fig. 2].

A large amount of plasma builds up due to the attraction of the magnetic field, and a current sheet is created by the magnetic field piling up. This corresponds to a strong increase in the plasma pressure which supports the current sheet and in the plasma temperature, as shown in the top and bottom right panels of Fig. 4. The non-adiabatic spike in the plasma temperature produced by the Ohmic heating $\eta J2$ (shown in Fig. 5) has an effect on the neutral temperature, which increases due to the thermal coupling between the species. As the neutrals are not completely coupled to the ions, they are expelled from the current sheet, as shown by the drop in the neutral density (bottom left panel of Fig. 4).

Both fluids are also heated up in the inflow region through frictional heating, the non-dimensional definition^{52} of which is

The frictional heating at *y* = 0 is shown in Fig. 5 compared to the Ohmic heating at *t* = 1.62 and *t* = 2.42. It increases with time during the first phases of reconnection, dropping only in the current sheet, where the Ohmic heating provides the major contribution by heating the plasma and increasing the plasma pressure and temperature.

### B. Current sheet and reconnection rate

Once the current sheet is generated, the MHD and PIP cases show very different reconnection processes. Laminar reconnection takes place in the MHD case, independent of the initial white noise perturbation without the onset of further instabilities. The long thin current sheet can be compared to the steady-state Sweet–Parker model. Both the length $\Delta MHD$ and the width $\delta MHD$ are estimated by taking the full width at 1/8 of the maximum current density *J _{z}*, respectively, along the

*y*-axis and the

*x*-axis. We choose this ratio as it accurately represents the termination of the reconnection region. The current sheet width and length are $\delta MHD\u223c0.02$ and $\Delta MHD\u223c1.49$. Using $\Delta MHD$ as characteristic length of the system and the maximum value of the Alfvén speed that occurs at the boundary of the current sheet, $vA\u223c3.26$, it is possible to calculate the Lundquist number. We find that $S=\Delta MHDvA/\eta \u223c9.7\xb7103$.

In the case of Sweet–Parker-like reconnection, the current sheet aspect ratio scales as $S\u22121/2$. From the value of *S* we obtain that the expected aspect ratio for the MHD case is $\delta /\Delta \u2009\u223c\u20091.01\xd710\u22122$, an estimate comparable to the measure obtained by $\delta MHD/\Delta MHD\u2009\u223c\u20091.57\xd710\u22122$. Having laminar reconnection in a long thin current sheet which does not develop instabilities nor break into smaller parts, we may suggest that the MHD case is subject to a reconnection process that is Sweet–Parker-like.

Unlike the MHD case, reconnection in the PIP case is nonlaminar. The presence of plasmoids breaks the current sheet into multiple thinner current sheets, which reconnect faster than the original structure. The dimensions of the current sheet at the time immediately before the generation of the first plasmoid (*t* = 2.42) are $\delta PIP\u223c0.01$ in width and $\Delta PIP\u223c0.46$ in length, estimated as $Jz,max/8$ along the *x*-axis and the *y*-axis, respectively. The onset of the plasmoid instability in a fully ionized plasma might take place below a critical aspect ratio^{29} of $1/200=5\xd710\u22123$: the same result was obtained for a multi-fluid plasma.^{42} Our current sheet aspect ratio is $\u223c1.9\xd710\u22122$, about four times larger than the predicted aspect ratio implying further physics may be involved in the onset of plasmoid formation. This is investigated in Sec. IV E.

In a partially ionized plasma, three different Alfvén speeds can be identified:^{17} a total Alfvén speed $vA,t$ related to the total density $\rho n+\rho p$, an ion Alfvén speed $vA,p$, based on the density of the charged particles, and an effective Alfvén speed $vA,e$, based on the combined density of charged particles and neutrals that are coupled through collisions. The expression for $vA,e$ might be non-trivial; however, it is possible to provide a close estimate for it from the plasma outflow velocity ($vout\u223cvA,e$). We chose to use the plasma velocity as the ionized fluid is the one directly accelerated by the reconnected magnetic field lines. The Alfvén speed is inversely proportional to the density: as *v _{A}* increases at the decrease in density, the ion Alfvén speed is bigger than the total Alfvén speed. The Lundquist number calculated by using $vA,t\u223c2.03$ is $St=vA,t\Delta PIP/\eta \u2009\u223c\u20091.75\xd7103$. If we consider the ion Alfvén speed only, which is $vA,p\u223c23.29$, the Lundquist number becomes $Sp\u2009\u223c\u20092.01\xd7104$, which is consistent with the threshold value $S=4\xd7104$ for the onset of the tearing instability and plasmoid formation.

^{21,27,42,67}

In the presence of collisional coupling, reconnection scales with the Lundquist number associated with the effective Alfvén speed $vA,e$. Estimating $vA,e\u223c9.09$ from the maximum outflow velocity, we find $Se\u2009\u223c\u20097.8\xd7103$, which is below the threshold value for the onset of the tearing instability. The discrepancy suggests that effects due to partial ionization might affect the dynamics of reconnection, allowing the formation of secondary plasmoids in the presence of a lower Lundquist number. The answer to this discrepancy between the models can be sought in the modifications due to the ion–neutral interaction.

The reconnection rate *M* is defined by

where $Jmax$ is the absolute maximum value of the current density inside the current sheet, *v _{A}* is the initial maximum bulk Alfvén speed ($vA,t$ in the PIP case), and $Bup$ is the initial maximum value of

*B*in the inflow. The MHD mean reconnection rate is $M=0.027\xb10.004$, and it appears to be higher than the Sweet–Parker rate by a factor of 2. The PIP reconnection rate, the average value of which is 0.08 ± 0.01, is approximately three times bigger than the MHD case and displays very sharp variations during the merging. Such fluctuations correspond to the formation and expulsion of the secondary plasmoids that disrupt the current sheet.

_{y}Figure 6 shows the temporal evolution of the maximum (black) and median (blue) drift velocity inside the current sheet during the reconnection phase at $t=[2.4,4]$ for the PIP reference case. Both are compared to the evolution of current density at the center of the current sheet ($x=0,y=0$), displayed in red. While the maximum drift velocity tends to oscillate quite drastically, specifically in connection of the major fluctuation in the current density, the median value tends to be approximately constant with a value of $\u223c0.1$, increasing smoothly in the last phases of the coalescence after *t* = 3.5, as the magnitude of *J _{z}* decreases. The peak in the median drift velocity is reached at the complete merging, where it reaches a value $|vD|\u223c1$.

The huge increase in the absolute maximum drift corresponds to the neutrals being expelled from the current sheet in the *x* direction during its collapse, and it takes place during the formation and expulsion of secondary plasmoids (for more details on secondary plasmoids see Sec. III D). This can be seen from the peak in the maximum drift velocity that is reached in the interval between the two central vertical dashed lines. These vertical lines represent the times identified in panels (b) and (d) of Fig. 2. The drop in the maximum $|vD|$ occurring approximately between *t* = 2.5 and *t* = 2.8 happens in correspondence of a relatively constant value of *J _{z}*. Such smooth variation of

*J*is linked to the formation and growth of the first secondary plasmoid ant the center of the current sheet. More details about the investigation of the current density are presented in Sec. IV B.

_{z}### C. Shocks

During reconnection and in the final phase of the merging in the MHD case [panels (c) and (e) of Fig. 2], there is evidence of shocks visible as thin elongated lines corresponding to both positive and negative peaks of the current density magnitude. The structures that can be distinguished in the current density are enhanced in the divergence of *v _{p}*, shown in Fig. 7, where they are identified as thin red lines.

The minimum in the divergence of the plasma velocity field identifies a region in which the flow is highly compressed, i.e., a shock. Across the shock, the magnetic field components *B _{y}* and

*B*drop, while plasma density and pressure rise steeply. The behavior of magnetic field and pressure identifies this as a slow-shock. The presence of slow-mode shocks is expected as they are a common feature in reconnecting systems, being part of a huge variety of fine structures that can be identified when plasmoid dynamics takes place.

_{z}^{68}

Comparing panels (c) and (d) of Fig. 2, the PIP case shows far fewer clear shock structures. Here we analyze the mechanisms suppressing slow-mode shocks in the PIP simulations. We examine the divergence of the plasma and neutral velocity fields (Fig. 8) at *t* = 3.32, the same time of panel (d) in Fig. 2. We focus on the divergence of the velocity field as large values of $\u2207\xb7v$ are a signature of shocks. Comparing it to the divergence of the velocity in the MHD case (Fig. 7), there are no structures that can be associated with slow-mode shocks. We present more information later on in Sec. IV B.

However, a wide range of structures appear in both $\u2207\xb7vn$ and $\u2207\xb7vp$, and they are particularly enhanced in the neutrals. In the plasma velocity divergence, the most prominent structure is associated with the neutral jet discussed in Sec. III E, but other structures cut the *x*-axis symmetrically at both sides of the reconnection region. These structures, which form in the neutrals and later couple to the plasma, are hydrodynamic shocks generated by the motion of the neutral species in the inflow and not slow-mode shocks as found in the MHD case.

During the merging, the neutrals are expelled and travel away from the reconnection region with a flow that is denser in the direction perpendicular to the current sheet. This can be seen from panel (a) of Fig. 14. In their motion, the neutrals interact with the dense plasma flow [panel (b) of Fig. 14], and they are halted by the collisions. The compression of the neutral flow leads to the formation of multiple shock fronts that are perpendicular to the *x*-axis. The neutral shocks coupling to the plasma manifest as the lines in $\u2207\xb7vp$. The lines, whose front moves away from the current sheet center, are present in the neutral *v _{x}* component, pressure and density but do not display a counterpart in the plasma variables.

Other hydrodynamic shocks are visible along the *y*-axis. These shocks are formed by the material accelerated inside a neutral jet, which will be examined in Sec. III E.

### D. Secondary plasmoids

During coalescence in the PIP case the central current sheet is subject to the tearing instability, and secondary plasmoids are produced as evident in panel (b) of Fig. 2. Figure 9 shows three secondary plasmoids forming, moving along the current sheet and being expelled. The motion along the current sheet is triggered by the white noise perturbation that breaks the symmetry of the system.

In the framework of plasmoid dynamics, it is interesting to investigate whether the secondary plasmoids have any characteristic in common with the initial plasmoids. We look at the force balance between the total pressure gradient and the Lorentz force ($J\xd7B\u2212\u2207p$), shown in Fig. 10 along the *y*-axis at *t* = 2.56 and *t* = 2.62. The vertical dashed lines are representative of the edges of the secondary plasmoids. At *t* = 2.56, plasmoid 1 can be identified at $y=[\u22120.1,0.15]$ and moves to the right at *t* = 2.62, while plasmoid 2 forms at $y=[\u22120.18,\u22120.1]$ at the later time. The force components cancel each other at the plasmoids location, while the current sheet around is still out of balance. Inside the plasmoids, the major contributions to the total force are provided by the magnetic pressure $B2/2$ (magenta) and the *y* component of the magnetic tension $(B\xb7\u2207)\xb7B$ (red), while the gradient of the neutral pressure (green) is negligible across the whole region. This suggests that the secondary plasmoids are in an almost force-free condition other than a small region at their center where the plasma pressure ($\u2212\u2207pp$ is shown in blue in Fig. 10) becomes significant.

The thermal coupling between ions and neutrals, whose terms are displayed on the right hand side of Eqs. (3) and (8), however, contributes to change the plasma pressure gradient with time. The effect of the thermal coupling is shown in Fig. 11, where the plasma pressure, the plasma temperature, and the two terms associated with the coupling with the neutrals (frictional heating and thermal damping) are displayed at the center of the secondary plasmoid 1 in the time interval $t=[2.56,2.64]$. The thermal damping term, whose non-dimensional definition is $(3/2\gamma )\alpha c(Tn,Tp,vD)\rho n\rho p(Tn\u2212Tp)$, drives the thermal equilibrium. When negative, the thermal damping indicates that energy is transferred from the hotter plasma to the neutrals. The plasma pressure decreases, together with the plasma temperature, under the effect of the thermal damping that reaches more negative values in time, a trend that is associated with energy passing from the plasma to the neutrals. The trend shown by the thermal damping reflects the neutral temperature, which tends to the plasma temperature until *t* = 2.62: after this time the plasmoid begins moving faster to the end of the current sheet as shown in Fig. 9. The frictional heating is also seen to increase remarkably with respect to time, as a result of the combined effect of the neutral pressure increasing from the inflow and the plasmoid motion. The coupling with the neutrals acts on $\u2207pp$ by sharpening its peak at the plasmoid center, but the gradient continues to be a relevant contribution to the total force. For this reason, the secondary plasmoids do not become completely force-free before they are expelled from the current sheet.

The detail of one of the secondary plasmoids reconnecting at one end of the current sheet is shown in Fig. 12, where *J _{z}*, the drift velocity magnitude, and the plasma and neutral velocity components are displayed. The plasma flow is faster than the neutrals, and this leads to a non-negligible drift velocity between the two fluids. Looking at both the current density and the drift velocity magnitude, a thin elongated vertical structure is observed from the current sheet to the center of the secondary plasmoid. This structure, visible in red in the bottom right panel of Fig. 12 (neutral

*v*), is a jet in the neutral flow, which extends in the direction opposite to the plasmoid motion, going back to the current sheet. This jet, accelerated by the neutral pressure gradient inside the plasmoid, is present only in correspondence of secondary plasmoids: there is no similar structure forming in the bigger coalescing plasmoids. The absence of this feature might depend on the fact that the bigger plasmoids are initially in a force-free condition and the pressure gradients are too small to expel the neutrals through a jet.

_{y}### E. Extended neutral reconnection jet

A prominent feature developing in the PIP simulation is the formation of a jet-like structure that extends asymmetrically along the *y*-axis during coalescence. This large structure must be distinguished from the small-scale neutral jets discussed in Sec. III D. In standard reconnection models, magnetic energy can be released to form a plasma jet, a feature that is also found in many observations.

From Fig. 13, however, we can clearly see that the neutral jet is significantly longer than the plasma jet. The ion velocity increases to supersonic values along the reconnection region, but the enhancement is localized near the center of the domain. The velocity of the extended neutral jet is supersonic, and the neutral Mach number reaches values of $\u223c1.6$. The larger drift velocity in Fig. 13 shows that the species are significantly decoupled in the jet.

A more detailed picture of the interaction of the two fluids along the jet structure can be provided by looking at the physical properties in the smaller region where the jet develops. Such region is identified in the top panel of Fig. 13. Plots of neutral and plasma densities and temperatures, frictional heating, and temperature difference between the two species are shown in Fig. 14.

The decoupling of neutrals and plasma along the jet is favored by the very low density of both species, as shown in panels (a) and (b) of Fig. 14. The two species reach similar peaks in temperature, as shown in panels (d) and (e) of Fig. 14, but the heating distribution is different for each fluid. The species are thermally decoupled, as shown by the difference between the neutral and ion temperatures in panel (f) of Fig. 14. During its evolution, the jet appears to be very turbulent. There is presence of many coherent vortices mostly concentrated at the jet truncation that are particularly evident in panels (c) and (f) of Fig. 14. Along the jet the ion temperature is the highest and reaches its maximum in correspondence of the current sheet, while the neutral temperature is higher than the ion temperature at the center of the vortices. Neutrals and ions are, however, heated up in the current sheet and along the jet by the thermal energy. The thermal energy is released through the frictional heating, defined in Eq. (23), which is associated with collisions between the two fluids, and it is shown in panel (c) of Fig. 14.

The velocity difference at the interfaces between the jet and the environment leads to the onset of shear flow instabilities. The sinusoidal shape of the jet is characteristic of the Kelvin–Helmholtz instability (KHI), a classical shear flow instability that tears apart vorticity sheets at the surface of separation of the two fluids.^{69–71} In order to confirm whether the system is KH unstable, we compare the neutral jet to the simple Bickley jet, a steady two-dimensional laminar jet which is unstable to the sinusoidal-mode of the KHI.^{69} Under the action of the KHI, the Bickley jet develops a sinusoidal structure at a preferred wavelength of $\u223c6.3$ times the characteristic flow half width.

The instability wavelength *λ* can be calculated at *t* = 3.16, where the instability is taking place from *y* = 0 and propagating downward along the jet. Measuring the wavelength as the distance between two vortices on the same side of the jet, $\lambda \u223c0.115$. The average half width of the jet, calculated as the distance between the peaks of maximum and minimum vorticity, is $\delta \u223c0.017$. We find an aspect ratio $\lambda /\delta $ of 6.6, which is similar to the value predicted for the Bickley jet undergoing KHI. We can conclude that the jet is subject to the KHI.

For our jet, the KHI is seen to evolve to a turbulent state. At the termination point, in this location shocks are generated and can be seen as weak structures in the neutral and drift velocities in Fig. 13 (top and bottom panels). These structures are the hydrodynamic shocks discussed in Sec. III C.

## IV. PARAMETER SURVEY

We investigate the changes in the coalescence process due to the diffusivity (*η*), the collisional coupling (*α _{c}*), the ion fraction (

*ξ*) and the plasma

_{p}*β*. In this section, we present a survey over these four key parameters of our physical system. The simulations are identified by numbers, and the respective physical parameters and the spatial resolution are listed in Table I.

Nr. . | Type . | Η
. | α
. _{c} | ξ
. _{p} | β
. | No. x grid points
. | No. y grid points
. | $\Delta x$ . | $\Delta y$ . |
---|---|---|---|---|---|---|---|---|---|

1 | MHD | 0.0005 | $\u221e$^{a} | 1^{a} | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

2 | PIP | 0.0005 | 100 | 0.01 | 0.1 | 6478 | 4862 | $1.2\xd710\u22123$ | $1.6\xd710\u22123$ |

3 | PIP | 0.0015 | 100 | 0.01 | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

4 | PIP | 0.005 | 100 | 0.01 | 0.1 | 1038 | 1550 | $3.9\xd710\u22123$ | $5.2\xd710\u22123$ |

5 | PIP | 0.015 | 100 | 0.01 | 0.1 | 1038 | 1550 | $3.9\xd710\u22123$ | $5.2\xd710\u22123$ |

6 | PIP | 0.05 | 100 | 0.01 | 0.1 | 1038 | 1550 | $3.9\xd710\u22123$ | $5.2\xd710\u22123$ |

7 | PIP | 0.15 | 100 | 0.01 | 0.1 | 1038 | 1550 | $3.9\xd710\u22123$ | $5.2\xd710\u22123$ |

8 | PIP | 0.5 | 100 | 0.01 | 0.1 | 1038 | 1550 | $3.9\xd710\u22123$ | $5.2\xd710\u22123$ |

9 | MHD | 0.0015 | 0^{a} | 0.01^{a} | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

10 | PIP | 0.0015 | 1 | 0.01 | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

11 | PIP | 0.0015 | 10 | 0.01 | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

12 | PIP | 0.0015 | 1000 | 0.01 | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

13 | PIP | 0.0015 | 3000 | 0.01 | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

14 | MHD | 0.0015 | $\u221e$^{a} | 1^{a} | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

15 | PIP | 0.0015 | 100 | 0.5 | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

16 | PIP | 0.0015 | 100 | 0.1 | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

17 | PIP | 0.0015 | 100 | 0.001 | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

18 | PIP | 0.0015 | 100 | 0.01 | 1 | 1038 | 1550 | $3.9\xd710\u22123$ | $5.2\xd710\u22123$ |

19 | PIP | 0.0015 | 100 | 0.01 | 0.01 | 1038 | 1550 | $3.9\xd710\u22123$ | $5.2\xd710\u22123$ |

Nr. . | Type . | Η
. | α
. _{c} | ξ
. _{p} | β
. | No. x grid points
. | No. y grid points
. | $\Delta x$ . | $\Delta y$ . |
---|---|---|---|---|---|---|---|---|---|

1 | MHD | 0.0005 | $\u221e$^{a} | 1^{a} | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

2 | PIP | 0.0005 | 100 | 0.01 | 0.1 | 6478 | 4862 | $1.2\xd710\u22123$ | $1.6\xd710\u22123$ |

3 | PIP | 0.0015 | 100 | 0.01 | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

4 | PIP | 0.005 | 100 | 0.01 | 0.1 | 1038 | 1550 | $3.9\xd710\u22123$ | $5.2\xd710\u22123$ |

5 | PIP | 0.015 | 100 | 0.01 | 0.1 | 1038 | 1550 | $3.9\xd710\u22123$ | $5.2\xd710\u22123$ |

6 | PIP | 0.05 | 100 | 0.01 | 0.1 | 1038 | 1550 | $3.9\xd710\u22123$ | $5.2\xd710\u22123$ |

7 | PIP | 0.15 | 100 | 0.01 | 0.1 | 1038 | 1550 | $3.9\xd710\u22123$ | $5.2\xd710\u22123$ |

8 | PIP | 0.5 | 100 | 0.01 | 0.1 | 1038 | 1550 | $3.9\xd710\u22123$ | $5.2\xd710\u22123$ |

9 | MHD | 0.0015 | 0^{a} | 0.01^{a} | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

10 | PIP | 0.0015 | 1 | 0.01 | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

11 | PIP | 0.0015 | 10 | 0.01 | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

12 | PIP | 0.0015 | 1000 | 0.01 | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

13 | PIP | 0.0015 | 3000 | 0.01 | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

14 | MHD | 0.0015 | $\u221e$^{a} | 1^{a} | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

15 | PIP | 0.0015 | 100 | 0.5 | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

16 | PIP | 0.0015 | 100 | 0.1 | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

17 | PIP | 0.0015 | 100 | 0.001 | 0.1 | 2062 | 3086 | $1.95\xd710\u22123$ | $2.6\xd710\u22123$ |

18 | PIP | 0.0015 | 100 | 0.01 | 1 | 1038 | 1550 | $3.9\xd710\u22123$ | $5.2\xd710\u22123$ |

19 | PIP | 0.0015 | 100 | 0.01 | 0.01 | 1038 | 1550 | $3.9\xd710\u22123$ | $5.2\xd710\u22123$ |

^{a}

These data are the effective values of the two-fluid parameters *α _{c}* and

*ξ*for the single-fluid cases, which are chosen as limits for the PIP simulations.

_{p}The PIP simulations in this section have the same resolution as the MHD cases or an even lower resolution ($\Delta x=3.9\xd710\u22123,\u2009\Delta y=5.2\xd710\u22123$). Due to the parameter variation changing the size of the central current sheet, it was possible to use a lower resolution without losing the possibility to resolve the current sheet.

### A. Variation of resistivity

We begin the investigation of partial ionization effects on the coalescence instability by considering the role of varying the resistivity in PIP simulations. The seven cases that are examined in this section ($\eta =$ 0.0005, 0.0015, 0.005, 0.015, 0.05, 0.15 and 0.5) are listed in Table I with the numbers from 2 to 8. The magnitude of *J _{z}* at the center of the current sheet (

*x*= 0,

*y*= 0) is displayed in Fig. 15, and it is seen to decrease as result of the increasing diffusion. The chosen location allows us to identify the beginning and the end of reconnection with the formation of the final plasmoid.

For $\eta =0.005$ (in black in Fig. 15) reconnection happens in the central current sheet with the formation of some secondary plasmoids, although fewer with respect to the less diffusive cases ($\eta =0.0005$, examined in Sec. III, and $\eta =0.0015$). The secondary plasmoids formation and expulsion are identified by a sudden drop (when the plasmoid is formed) followed by a drastic increase back to previous values of $|Jz|$ as soon as it moves away from the center of the current sheet. In the simulation with $\eta =0.015$ the current is highly diffused. During the plasmoid merger there is no observed formation of secondary plasmoids. This is the first case in which the diffusion is high enough to prevent the formation of secondary plasmoids during coalescence in a PIP case. A similar result is obtained for the case with $\eta =0.05$. The cases having the highest resistivity ($\eta =0.15,0.5$) do not develop the coalescence instability, as the initial plasmoids are quickly diffused. As *η* is increased, the Lundquist number *S* decreases: when *S* becomes sufficiently low, the tearing instability stops taking place in the central current sheet, thus explaining the lack of secondary plasmoids.

The variation of resistivity also plays a role in changing the time scale of coalescence. For higher *η*, coalescence starts at a later time in each simulation. This is shown in Fig. 15 by the position in time of the first drops in *J _{z}*. Such effect depends on how efficiently the current sheet is generated in the phase where the two plasmoids are approaching, as for smaller

*η*it is easier to buildup current.

### B. Variation of collisional frequency

Here we investigate the effects of the collisional coupling between ions and neutrals. We compare the simulations listed in Table I with the numbers 3, 9, 10, 11, 12, 13, and 14. An MHD case with total density equal to $\rho p=1$ (simulation 14) is taken as the limit case $\alpha c\u2192\u221e$. A second MHD case (simulation 9) having initial density equal to $\rho p=0.01$ is considered as the limit for $\alpha c=0$ (i.e., no collisions), as this is the same plasma density set for the PIP cases. For $\alpha c=\u221e$ the neutral and plasma species are entirely coupled hence the system behaves like a single fluid MHD model, with the density and pressure being the bulk (ion + neutral) values. For $\alpha c=0$, the species are completely decoupled, i.e., the plasma evolves independently from the neutrals and can be considered to be a single fluid MHD system with the density/pressure based on the plasma values only. An equivalent MHD simulation can be performed by changing the initial plasma beta. While the magnetic field strength is unchanged, the variation of plasma beta modifies the pressure. To maintain the same initial temperature of the calculation, we use a lower density (the same as the plasma density of our two-fluid calculations). Therefore, the difference in the plasma density between the two limit MHD cases result in an effective difference in the plasma *β* of the two simulations.

In Fig. 16, we display the variation of current density *J _{z}* at $x,y=0$ with respect to time. The beginning of reconnection is identified with the first minimum occurring in the current density, when the current sheet is compressed the most by the two plasmoids. Decreasing

*α*the time scale for the plasmoid coalescence becomes considerably shorter. The time scale does not vary linearly, but it shows the presence of two accumulation points, corresponding to the two limits identified by the MHD simulations. At lower

_{c}*α*($\alpha c=1,10$), the time scale tends to approach the one for $\alpha c\u21920$, while at higher

_{c}*α*($\alpha c=1000,3000$) coalescence takes place in a time interval similar to the one obtained for $\alpha c\u2192\u221e$.

_{c}The more ions and neutrals are coupled, the later reconnection starts. The slowing down of the first phase prior to the onset of reconnection can be associated with the damping effect that neutrals have on the ions in the inflow, which increases at the increase in *α _{c}* as the ions (which fraction is much smaller than the neutral fraction) interact with a higher number of neutrals. This result is consistent with previous simulations.

^{50}We can understand this result by comparing the ion Alfvén time $\tau A,p$ with both the ion and the neutral collisional times. This leads to the identification of two values for

*α*that define when the species couple with each other through collisions. For $\alpha c\u223c20$, the plasma–neutral collisional time $\tau col,pn=(\alpha c\rho n)\u22121$ becomes smaller than $\tau A,p$ and the ions couple with the neutrals. The coupling of the neutrals on the ions, that happens when the neutral–plasma collisional time $(\alpha c\rho p)\u22121$ becomes smaller than $\tau A,p$, takes place for $\alpha c\u223c2000$.

_{c}As *α _{c}* increases, reconnection takes a longer time before the plasmoids are completely merged. This results in a decrease in the average reconnection rate as

*α*increases, as shown in the top panel of Fig. 17, where the error bar is given by the standard deviation. The reduced reconnection rate can be explained by the variation of the effective Alfvén speed, as introduced for the two-fluid case in Sec. III B. If the collisional frequency is increased, the effective density of the magnetic fluid increases as more neutrals interact with the ions and consequently $vA,e$ decreases. The variation of the Lundquist number is shown in the central panel of Fig. 17, as calculated from $vA,e$.

_{c}After reconnection begins, the PIP cases (Fig. 16) show a strong fluctuation of the current toward less negative values, behavior, that is, not present for a full coupling (MHD case for $\alpha c\u2192\u221e$). This fluctuation is a consequence of plasmoid formation in the current sheet. The growth and expulsion of plasmoids contributes initially to slow down the reconnection by saturating the negative current in a blob between the merging plasmoids and then to leave the current sheet unstable, leading to the formation of further plasmoids. For the simulations with secondary plasmoids, *S* is calculated at the time before the formation of the first plasmoid in the central current sheet, while in the case of simulations without secondary plasmoid *S* is evaluated at the first minimum of *J _{z}* at the center of the current sheet. The value of the Lundquist number decreases at the increase in the collisional coupling, following the variation of both the effective Alfvén speed and the characteristic length of the current sheet.

*S*is more sensitive to the variation of $vA,e$, as such factor varies over an interval covering orders of magnitude, while Δ, being about the same size of the plasmoids, displays a lesser variation in length across the cases. As shown in Fig. 17,

*S*appears to be above the threshold of 10

^{4}for all the cases showing formation of secondary plasmoids, with the exception of the case with $\alpha c=1000$ and 3000. The effects that might lead to the sub-critical plasmoid formation we have found here will be discussed in Sec. IV E.

Looking at the results of Sec. III C, slow-mode shocks are apparently absent in the PIP case with $\alpha c=100$. However, slow-mode shocks can be generate in two-fluid environments.^{72} Therefore, we want to investigate whether slow-mode shocks are produced at a higher or lower collisional coupling that approach the MHD cases, and whether they form but are dissipated at a later time of coalescence. Slow-mode shocks are indeed generated in the two-fluid simulations, and their propagation is clearly visible particularly at low *α _{c}*, as shown in the first two panels (simulations 10 and 11) of Fig. 18. At the increase in collisional coupling ($\alpha c\u226510$), they are damped and disappear as a consequence of two-fluid effects. At higher collisional frequencies ($\alpha c=100,1000$ displayed in the last two panels of Fig. 18), the turbulent motion set by the neutral jet and the propagation of hydrodynamic shocks disrupt the slow-mode shock front in proximity of the reconnection region. The presence of hydrodynamic shocks, absent for $\alpha c\u22641$, increases with the coupling between the two species around the inflow region as the two species that move in opposite directions are subject to an increasing interaction. For $\alpha c=100$, the slow-mode shocks cannot be detected. Due to the better coupling, in the case with initial $\alpha c=1000$, the hydrodynamic shocks show a similar behavior as the slow-mode shocks: their front moves along the magnetic field lines, as shown by their position with respect to the separatrices in Fig. 18.

### C. Variation of ion fraction

In this section, we compare five cases with a different ion fraction *ξ _{p}* in the range $10\u22123\u22125\xd710\u22121$, corresponding to the numbers 3, 14, 15, 16, and 17 listed in Table I.

The variation in time of *J _{z}* at the center of the current sheet is displayed in Fig. 19, with the speed in the process of coalescence drastically increased at the decrease in

*ξ*. Such behavior, which shows a variation in the time scale of both ideal phase (when the plasmoids attract each other and the current sheet is formed) and reconnection phase, might be explained by similar arguments to those presented in Sec. IV B.

_{p}At the variation of the ion fraction, the initial ion Alfvén time scale $\tau A,p$ and the ion collisional time scale increase with *ξ _{p}*, while the neutral collisional time scales decreases. We compare $\tau A,p$ with the ion and the neutral collisional times to find the values of

*ξ*for which the two species couple with each other. In our range, the ion collisional time $(\alpha c\rho n)\u22121$ is always smaller than $\tau A,p$ with the only exception of the MHD case (simulation 14) which is taken into account as the limit value for $\xi p=1$. Therefore, the ions are always coupled to the neutrals in all the PIP simulations in this section. The value of

_{p}*ξ*below which the ion dynamics becomes fast enough to decouple from the neutrals is $\u223c4\xd710\u22124$. Conversely, the neutrals coupling on the ions takes place for $\xi p=0.075$. In the PIP cases with the highest ion fractions ($\xi p=0.1,0.5$), in which the magnetic forces are felt by a significant portion of the fluid, the coalescence develops in a similar way as in the MHD case.

_{p}The reconnection rate, shown in the top panel of Fig. 20, decreases as the ion fraction increases, following the variation in the coalescence time scale. At the increase in *ξ _{p}* the effective Alfvén speed $vA,e$ decreases, following the increase in the plasma density. This affects the Lundquist number (central panel of Fig. 20), which in turn extends the time scale of the reconnection phase. Sub-critical plasmoid formation, discussed in Sec. IV E, is observed for the case at the lowest ion fraction ($\xi p=0.001$), with

*S*below the threshold limit of 10

^{4}at the onset of the tearing instability.

### D. Variation of plasma beta

In this subsection, three cases at different plasma *β* (simulations 3, 18 and 19 in Table I) are investigated. The speed of coalescence is greatly affected by the change of *β*, and increases at the decrease in this parameter as shown in Fig. 21. This is explained by the fact that at smaller *β* the magnetic field get stronger with respect to the plasma pressure: the Alfvén speed increases (as its numerator increases), and so does the Lundquist number.

The consequences of plasma *β* variations are not directly associated with a two-fluid effect, as the different dynamics depends uniquely on the variation of **B** and can be purely reduced to an MHD effect. However, it provides important context for how changing different parameters in the two-fluid case, which effectively alter *v _{A}*, can increase the merger rate.

The presence of sub-critical plasmoid formation is observed for the case with the largest plasma *β* (*β* = 1), where the formation of a single secondary plasmoid takes place for a Lundquist number smaller than 10^{4}. The physics behind this case is discussed in Sec. IV E.

### E. Sub-critical plasmoid formation

Several PIP simulations presented in Secs. IV A–IV D (identified in Table I by number 2, 4, 12, 13, 17, and 18) show the signs of secondary plasmoids formation in the central current sheet below the Lundquist number threshold $S=104$. The values of the critical *S* found in these cases are consistent with the results obtained by recent studies of magnetic reconnection in a multi-fluid partially ionized plasma,^{38} already discussed in Sec. I. One might expect that the lower critical Lundquist number in our simulations could depend on the initial setup, and, in particular, the low ion fraction. In terms of the initial setup, a comparison with previous works is difficult: many studies^{42,43} are performed with a static current sheet setup, while our simulations present a driven reconnection, which leads to a different current sheet structure and dynamics. A more detailed evaluation of the role played by the initial perturbation and the random white noise in Eq. (22) on the onset of the tearing instability must be investigated in future developments of this research. Several studies already found that the role of the amplitude of perturbation noises^{37,39–41} is major in determining the critical Lundquist number and current sheet aspect ratio in a range of initial configurations. However, tests performed on the variation of the white noise perturbations on our simulations proved that the secondary plasmoids are always generated at the same time in every simulation that has the same initial set of parameters.

In this section, we often refer to the effective Lundquist number as the critical Lundquist number for simulations displaying secondary plasmoid formation. The critical Lundquist number for the onset of the tearing instability was not isolated in previous studies,^{29,38} but it was considered to fall in an interval of values whose limits are set by the absence (lower boundary) and presence (upper boundary) of plasmoid formation. Such intervals had been determined by varying two main parameters of the current sheet: characteristic length and resistivity. In our study, this interval is determined by changing the characteristic length scale of the current sheet. We do not change *L* directly, but it varies in time as the system evolves. The time interval between two outputs is small enough that the current sheet length does not display a large variation: this reduces to have Lundquist numbers that are very close to each other before and after the first secondary plasmoid appear, which we identify by the formation of an *O*-point in the current sheet magnetic field. For such reason, we choose to identify the critical Lundquist number as a single value rather than specifying an interval.

We might expect the sub-critical plasmoid formation to be triggered in all the PIP cases because of the small ion fraction. However, the critical Lundquist number obtained for the MHD case at lower ion density (matching that of the plasma density of the PIP simulations), shown in Fig. 17, is well above the threshold of 10^{4} (limit case for $\alpha c=0$, listed with number 9 in Table I). This proves that the change in the plasma density itself, along with the variation of the plasma *β* of the plasma that this change implies, does not affect the system so that it develops sub-critical plasmoids. We also see that for the MHD calculation, that is, equivalent to $\alpha c=\u221e$ the current sheet is stable for $S=9.7\xd7103$, again implying a critical $S>104$. At the same time, PIP cases characterized by an *S* and an aspect ratio $\delta /\Delta $ incredibly similar to this MHD calculation (see cases 12 and 13, where *α _{c}* is, respectively, 10

^{3}and $3\xd7103$, in Fig. 17), show sub-critical plasmoid formation. The inclusion of the coupling between the ion and neutral fluids (effectively looking at systems between those two MHD simulations) allows plasmoids to form below $S=104$, even without changing the plasma

*β*, the initial conditions or the perturbation magnitude. Therefore, such sub-critical plasmoid formation does not depend on the characteristics of the plasma itself, but it might be ascribed to the combined result of two-fluid effects.

From previous studies of the onset of the tearing instability in partially ionized plasmas,^{73} the critical aspect ratio for the initiation of a growth rate independent of the Lundquist number and neutral to ionic density was derived for three regimes (coupled by collisions, intermediate and uncoupled). In the intermediate regime, where ions and neutral are partly coupled, the critical aspect ratio for a generic equilibrium configuration is

where the collisional frequency *ν _{pn}* is equal to $\alpha c\rho n$ and $\tau A,p=\Delta /vA,p$ is the ion Alfvén time scale. In general, this critical aspect ratio is bigger than both the threshold $S\u22121/3$ obtained from classical arguments that involve the Lundquist number and the value of 1/200 found in several works.

^{29,42}

Comparing our current sheet aspect ratios at the time of the onset of the tearing instability in the PIP simulations with the critical aspect ratio, we find that in most of the cases the aspect ratio is still smaller than the critical value obtained by Eq. (25). Therefore, the sub-critical plasmoid formation cannot be explained uniquely by applying such condition for the onset of the tearing instability.

Focusing on a physical explanation for the sub-critical plasmoid formation, we take a close look at what happens in the current sheet before the onset of tearing instability. We examine in detail the PIP case corresponding to simulation 13 in Table I and evaluate how the two-fluid parameters vary up to the time immediately before the formation of a central plasmoid, which signature appears in the current density and magnetic field at *t* = 4.7. This PIP case is characterized by an initial strong coupling between the two species as the initial $\alpha c(0)$ is set equal to 3000, and the Lundquist number before the onset of the tearing instability is $S=3\xd7103$.

The plasma velocity along *x* is shown in row (a) of Fig. 22. As enhanced by the contour lines, before the formation of the secondary plasmoid the plasma in the inflow slows down around the center of the current sheet where a stagnation point in the flow exists, forming a pinching flow that promotes reconnection in two points above and below the current sheet central point.

In order to understand how the plasma motion is slowed down, we evaluate the force balance in the inflow region. All the force contributions are in equilibrium with each other, with the exception of the drift force, stated on the right hand side of Eqs. (2) and (7), *x* component of which for the plasma is shown in row (b) of Fig. 22. The drift force magnitude appears to be quite large; however, this does not mean that the drift velocity between the two species is also big in the inflow region. The drift force has a strong dependency on the collisional term *α _{c}*, which is itself dependent on both the drift velocity and the temperatures of the two fluids. Because of the strong dependence on the

*α*term, when the species are weakly coupled (and consequently the drift velocities are large), the drift force is small.

_{c}The drift force acts by slowing down the flow as it approaches the current sheet: as its sign depends on the term $(vn\u2212vp)$, its direction suggests that the plasma is slowed down by the interaction with the neutrals, whose velocity in the inflow is slower than the plasma. The neutrals are naturally slower as they move toward the center of the current sheet as they are being pulled in by the plasma motions. However, as shown in panel (c), the neutral pressure increases at the center slowing the inflow motions and increasing the drag. The neutral pressure increases under the effect of the Ohmic heating $\eta J2$ and of the adiabatic compression of the neutrals in the inflow. As already discussed for the reference PIP case in Sec. III A, the Ohmic heating plays a major role in heating the plasma inside the current sheet (see Fig. 5), which results in the increase in the plasma pressure and temperature. As the two species are thermally coupled, the Ohmic heating acts by indirectly heating the neutrals, driving their pressure to increase. The non-adiabatic contribution from the Ohmic heating is responsible for the increase in the plasma temperature inside the current sheet, which affect the neutral temperature. The heated neutrals expand in the inflow region, leading to an adiabatic compression directed out of the current sheet, which contributes to increase the neutral pressure.

The sub-critical plasmoid formation can therefore be triggered by the two-fluid interaction between plasma and neutrals. The PIP cases at lower collisional coupling and higher ion fraction, however, develop the formation of critical plasmoids following the linear condition for the tearing instability before this non-linear mechanism becomes important. Therefore, the two-fluid interaction plays a lesser role in the formation of plasmoids when the effect of neutrals is weak, i.e., at smaller *α _{c}* and higher

*ξ*.

_{p}## V. DISCUSSION

Among the processes promoting the development of fast magnetic reconnection, the coalescence instability can play an important role by driving the interaction of plasmoids (and their subsequent reconnection) on dynamic time scales. We have investigated how plasmoid coalescence behaves in a partially ionized plasma, a situation reflected in a range of solar atmospheric layers and, in particular, the solar chromosphere. Through the comparison of a fully ionized case (MHD) and a partially ionized case (PIP), we find that partial ionization noticeably shortens the coalescence time scale and creates new dynamics, producing neutral jets and secondary plasmoids, suppressing slow-mode shocks while promoting hydrodynamic shocks, and leading to sub-critical plasmoid formation.

Observations of chromospheric anemone jets performed with the Ca II H filtergram of the Solar Optical Telescope onboard Hinode^{8} showed the presence of recurrent plasmoids expulsion with a size of about a few hundred km from the jets. Assuming an approximate diameter $\u2205MAX\u223c500$ km our characteristic length is $\u223c100$ km, as in our system the initial plasmoid length is equal to 4*L*. Knowing that the sound speed is about 10 km s^{−1}, we identify a time scale of $\u223c10$ s. Taking the appropriate ion–neutral collisional cross section for elastic scattering and charge exchange^{74} and the characteristic temperature, density and ion fraction^{18,45–47} of the chromosphere the ion–neutral collisional frequency^{75,76} is between 10^{3} and 10^{6} s^{−1}. In non-dimensional form, to compare with our simulations, an equivalent *α _{c}* to the one observed for the Hinode plasmoids is in the range of $104\u2212107$.

A prediction on the behavior at these higher *α _{c}* can be made by studying the simulations presented in Sec. IV B, whose trend is shown in Figs. 16 and 17. When

*α*increases, the evolution tends to the one of the upper limiting MHD case ($\alpha c\u2192\u221e$): such case can be considered as an accumulation point. In the evolution of the current density, all the simulations at the chromospheric

_{c}*α*would reach the beginning of reconnection in a time interval between the current minimum of the PIP case with $\alpha c=3\xd7103$ and the one of the MHD case, so it is sensible to say that the coalescence evolution and time scale can be predicted. These values suggest that the coalescence between the biggest observed plasmoids would take place in a regime, that is, almost MHD. At the specific values for ion fraction, plasma

_{c}*β*, and resistivity used in this study, faster coalescence would become relevant for plasmoids with a diameter of $\u223c3$ km. Such length scale can be found by comparing the highest collisional frequency of our simulations with the dimensional chromospheric collisional frequency. Taking 1 km as a characteristic length and a sound speed of 10 km s

^{−1}, we obtain a time scale $\tau =0.1\u2009s$. This leads to a collisional frequency of $\u223c104$ s

^{−1}, consistent with the collisional frequency observed in the chromosphere and with our case at a non-dimensional $\alpha c=3\xd7103$.

This length scale is, however, dependent on the parameters of the medium studied. In many regions of the solar chromosphere the resistivity is often smaller, which results in the formation of more plasmoid dynamics, and both *ξ _{p}* (equal to 0.01 in the

*α*section of the parameter survey) and the plasma

_{c}*β*(set equal to 0.1) can be lower in the chromosphere, leading to an enhancement of two-fluid effects. At lower

*ξ*and

_{p}*β*, for the range of

*α*considered in this study, PIP effects that include a faster coalescence and the sub-critical plasmoid formation would become important for larger plasmoids and are potentially observable on scales that are currently resolved. The 3 km plasmoid length scale must therefore be considered as a lower limit for PIP effects to become observable.

_{c}Partial ionization plays a role in changing the way many effects, such as the Hall diffusion, develop and act on magnetic reconnection. The Hall effect results from the different velocities of electrons and ions,^{77} is dependent on the ion fraction and is seen to play a major role in both highly ionized and weakly ionized plasmas.^{78} In a partially ionized environment, the physical scales over which the Hall diffusion becomes important drastically change from the ones found for fully ionized plasmas.^{78} Several works^{46,47,77} show that in magnetic flux tubes, of which plasmoids are to be considered 2D sections, the Hall diffusion is less important than the contribution from ion–neutral effects at the chromospheric heights. As such we viewed that including the Hall effect is beyond the scope of this paper.

A very important effect in evidence in our results is the sub-critical plasmoid formation discussed in Sec. IV E. Many studies on the onset of the tearing instability in partially ionized plasmas often focused on linear stability criteria,^{33,34,36,48,73} while in our study plasmoid formation is also promoted by a nonlinear effect linked to the two-fluid collisional and thermal coupling. The role played by the coupling between ions and neutrals in determining the dynamics of plasmoid formation has already been acknowledged in a previous study^{48} looking at the onset of the tearing mode at several levels down to the kinetic scale. In that study the authors only use a linear criterion for the onset of the tearing instability, while in this work we find that nonlinear effects are able to play a major role in plasmoid formation. Including the nonlinear physics may result in new physical parameter regimes that are able to tear down to kinetic scales.

The physics that leads to sub-critical plasmoid formation in our simulations is expected to be largely affected by the non-equilibrium ionization–recombination processes. These are currently not included in our model as the interaction between the two species is provided uniquely by elastic collisions and charge-exchange, but future developments of this research strictly require their inclusion. As pointed out from previous studies of magnetic reconnection in a multi-fluid partially ionized plasma at low plasma *β,*^{38,42,43,79–81} the non-equilibrium ionization–recombination effect leads to a strong ionization of the material in the reconnection region. In a recent paper^{80} magnetic reconnection has been examined through 2.5D simulations in weakly ionized plasmas with an initial *ξ _{p}* ($\xi p=10\u22124$) lower than our reference cases. Their results suggest that, for a plasma

*β*larger than 1 and weak reconnection magnetic fields, the non-equilibrium ionization–recombination effect is responsible for a strong ionization of the neutral fluid in the reconnection region and a faster reconnection rate occurring before the onset of plasmoid instabilities. These are consistent with a previous paper,

^{43}where an increase by an order of magnitude was recorded for the ionization degree within the current sheet during the reconnection process. The strong ionization is responsible for a bigger interaction between the neutral fluid and the plasma, which will be better coupled both in the inflow and outflow region. These very same effects are to be expected in the case discussed in Sec. IV E: the drastic increase in temperature due to the Ohmic heating in the reconnection region would promote the ionization of the neutral fluid that in the absence of such process is forced to expand outward, halting the plasma inflow. In the case of a small plasma

*β*smaller than 1 like in our study, however, plasmoid instability is still the main process promoting fast magnetic reconnection.

^{80}

The ionization/recombination effects are largely affected by the action of the ionization potential. When collisional ionization takes place, the work done against the ionization potential to ionize the atom removes energy from the electrons and acts as a cooling term in the plasma.^{82} As the recombination process is associated with photons being released, this overall effect can be modeled as a radiative loss. Researches investigating the role of radiative cooling in magnetic reconnection^{83,84} proved that such process, linked to the collisional ionization, thins the reconnection layer by decreasing the plasma pressure and density inside the current sheet. Therefore, the inclusion of radiative losses speeds up reconnection to rates that are bigger than the ones found in models without radiation and might lead to time scales and outflows that are consistent with those found in spicules and chromospheric jets.^{66} In a recent study,^{66} it was found that a strong recombination process in the reconnection region, combined with Alfvénic outflows, can lead to a fast reconnection rate independent of Lundquist number. While the decoupling of neutrals and plasma has been recorded in the inflow region, these findings show that the two fluids are well coupled in the outflows, which is opposite to what we find for our intermediate coupled case (see Sec. III E). We, therefore, expect that including the radiative losses by adding an ionization potential would lead to a better coupling of the two species around the reconnection region as well as lower temperatures and a faster reconnection rate.

Radiation is not only important in terms of the radiative losses from the ionization/recombination processes in the chromosphere. Beyond the action of the ionization potential, effects of ionization and excitation could be triggered by an external radiation field. This is a fundamental factor in the chromosphere, where the ionization degree is largely determined by the incident external radiation.^{85–88}

In a recent study^{38} it was demonstrated that the plasmoid cascading process, for which the current sheet breaks into smaller sections following the formation of multiple plasmoids, is terminated in the MHD scale. The progressive reduction of secondary plasmoids in the PIP cases as an MHD-like regime is approached is an aspect that has been already marginally observed in our simulations, specifically in the parameter survey linked to the variation of the collisional coupling (Sec. IV B). Multiple plasmoids are seen to form in the simulations having an initial lower *α _{c}*, while at higher collisional coupling only two and one secondary plasmoids are produced, respectively, for $\alpha c=1000$ and 3000. Moreover, as discussed in Sec. IV E, the secondary plasmoids formed at higher

*α*are generated by the pinching action of the neutrals in the inflow, rather than the onset of instabilities in the current sheet. Therefore, we expect the plasmoid cascading process to be interrupted in these cases approaching the MHD regime. A more detailed study needs to be performed on these secondary plasmoids number and characteristics, however, we can already confirm a good agreement with the trend showed by previous studies.

_{c}^{38}

## ACKNOWLEDGMENTS

The authors are grateful to Dr. N. Nakamura for the inspiration of his original work on this problem that led to this study. A.H. and B.S. are supported by STFC Research Grant No. ST/R000891/1. A.H. is also supported by his STFC Ernest Rutherford Fellowship Grant No. ST/L00397X/1.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request. The (PIP) code is available at the following url: https://github.com/AstroSnow/PIP.

## References

*Proceedings from International Astronomical Union Symposium No. 6*, edited by