Collisional ionisation and recombination effects on coalescence instability in chromospheric partially ionised plasmas

Plasmoid-mediated fast magnetic reconnection plays a fundamental role in driving explosive dynamics and heating, but relatively little is known about how it develops in partially ionised plasmas (PIP) of the solar chromosphere. Partial ionisation might largely alter the dynamics of the coalescence instability, which promotes fast reconnection and forms a turbulent reconnecting current sheet through plasmoid interaction, but it is still unclear to what extent PIP effects influence this process. We investigate the role of collisional ionisation and recombination in the development of plasmoid coalescence in PIP through 2.5D simulations of a two-fluid model. The aim is to understand whether these two-fluid coupling processes play a role in accelerating reconnection. We find that in general ionisation-recombination process slow down the coalescence. Unlike the previous models in G. Murtas, A. Hillier \&B. Snow, Physics of Plasmas 28, 032901 (2021) that included thermal collisions only, ionisation and recombination stabilise current sheets and suppress non-linear dynamics, with turbulent reconnection occurring in limited cases: bursts of ionisation lead to the formation of thicker current sheets, even when radiative losses are included to cool the system. Therefore, the coalescence time scale is very sensitive to ionisation-recombination processes. However, reconnection in PIP is still faster than in a fully ionised plasma environment having the same bulk density: the PIP reconnection rate ($M_{_{\operatorname{IRIP}}} = 0.057$) increases by a factor of $\sim 1.2$ with respect to the MHD reconnection rate ($M_{_{\operatorname{MHD}}} = 0.047$).


I. INTRODUCTION
Magnetic reconnection takes place in a wide range of astrophysical settings 1 : it occurs in presence of a non negligible resistivity, when magnetic field lines change their connectivity altering the topology of the magnetic field. During reconnection, the frozen-in constraint imposed by ideal magnetohydrodynamics (MHD) no longer applies and the energy released by reconnection energises particles to high energies and heats the plasma 2,3 . Explosive events in the solar chromosphere, such as chromospheric jets 4,5 and Ellerman bombs 6,7 , are believed to be driven by fast magnetic reconnection, responsible of efficiently releasing the stored magnetic energy into thermal and kinetic energy 2,3 at shorter time scales than the classical reconnection models.
The hypothesis of fast magnetic reconnection as a driver of explosive phenomena in chromospheric plasmas has been reinforced by the identification of plasma blobs in the outflow of chromospheric jets and UV bursts [8][9][10] , colliding and merging with each other before being ejected from the current sheet. These structures are generally interpreted to be plasmoids, concentrations of current density trapped in closed loops of magnetic field lines that are commonly present in reconnecting systems [11][12][13][14][15][16][17] . Plasmoids are believed to play a major role in speeding up reconnection to the time scales found in observations, as their formation directly affects the current sheet size 18 . Plasmoids can be the result of the tearing instability, which breaks thin current sheets (where the current sheet aspect ratio δ /L 1) into fragments 11,15,19 . The resulting high current densities in each of these fragments facilitate a high reconnection rate 20 . a) Electronic mail: gm442@exeter.ac.uk Plasmoid formation due to the instability of current sheets has been extensively examined through numerical studies 17,[21][22][23][24][25][26][27] . This mechanism is dependent on the value of the Lundquist number 13,28,29 , S = Lv A η , where L is a characteristic length of the system, v A is the Alfvén speed and η is the diffusivity. Several works proved that it is possible to speed up reconnection in thin current sheets for a critical value of the Lundquist number typically ∼ 10 3 − 10 4 for fully ionised plasmas 14,16,17,20,27,[30][31][32][33][34][35] . Above this limit, current sheets become unstable and plasmoids formation occurs.
Once plasmoids are generated, they are pulled against each other and merge, in a mechanism that further increase the reconnection rate: this process is called coalescence instability 2,36 . Plasmoid coalescence occurs in the nonlinear tearing mode phase between plasmoids sharing an X-point and evolves through two separate phases. The first is an ideal MHD phase with a growth rate that is almost independent of η 36 . In this phase, the two plasmoids move close and a current sheet forms between them. The second part of coalescence evolves in a resistive phase, where the current sheet begins to reconnect 2 . As further plasmoids are formed in between the coalescing plasmoids, this instability can become a fractal process operating at multi-spatial scale, depending on the Lundquist number 13,28,29 .
A large portion of plasmas across the universe are only partially ionised, for example the chromospheric plasma, whose ionisation degree falls in the range 10 −4 − 10 −1 , as reported by many studies 20,[37][38][39][40] . The presence of neutral species might alter the plasmoid dynamics, as further physical processes develop from the coupling with charged particles. The role of partial ionisation on the onset of magnetic reconnection and resistive tearing instability was investigated in several studies 8,18,20,[41][42][43][44] . It has been observed that partial ionisation largely modifies the reconnection rate by changing the Alfvén speed, which in turn affects the Lundquist number of the system 42 and the conditions for plasmoids formation.
The physics of a reconnecting current sheet is also affected by processes of collisional ionisation 45 and recombination, which actively modify the relative abundance of the plasma species. As pointed out from many studies of magnetic reconnection in a multi-fluid partially ionised plasma at low β 7,42,43,46-48 the non-equilibrium ionisationrecombination leads to a strong ionisation of the material in the reconnection region and a faster reconnection rate developing before the onset of plasmoid instabilities 48 . Past studies 43 reported an increase of the ionisation degree by an order of magnitude within the current sheet during reconnection. The strong ionisation is responsible for a larger interaction between the neutral fluid and the plasma, with a stronger coupling occurring both in the inflow and outflow region. In case of a plasma β smaller than 1, however, plasmoid instability remains the main process promoting fast magnetic reconnection 48 .
The ionisation-recombination effects are further enhanced by the action of the ionisation potential. When collisional ionisation takes place, the energy expended by a free electron to release a bound electron results in a net loss of energy from the plasma 49 . As the recombination process is associated with changes in energy levels and photons being released, this overall effect can be modeled as a radiative loss. Studies investigating the role of radiative cooling in magnetic reconnection 50,51 proved that the inclusion of the ionisation potential thins the reconnection layer by decreasing the plasma pressure and density inside the current sheet. Therefore, adding radiative losses speeds up reconnection to higher rates than the ones of models without radiation and might lead to time scales and outflows that are consistent with those found in spicules and chromospheric jets 52 .
Relatively little is known about how the coalescence instability is altered by the multi-fluid physics of a partially ionised plasma. In our preliminary study 53 we found that, in an ionneutral plasma with fluids coupled through elastic collisions, partial ionisation speeds up both phases of coalescence and promotes non-linear effects during reconnection. In this paper we examine in detail the effects of collisional ionisation, recombination and optically thin radiative losses on the coalescence instability developing in partially ionised plasmas. This model improves upon our previous work 53 by including the contribution of these processes. We aim to understand how the new type of coupling between charges and neutral species influence the reconnection rate, compared to our previous research. In Section II we discuss the two-fluid model employed for our simulations. In Section III we report the results of our 2.5D simulations. The results are discussed in Section IV.

II. METHODS
Numerical simulations are performed of the coalescence instability in 2.5D using the (PIP) code 54 . The process is studied in single-fluid fully ionised cases (MHD) and twofluid partially ionised cases (PIP). MHD cases are modelled by a charge-neutral ion-electron hydrogen plasma, while PIP cases are characterised by a neutral fluid and a hydrogen plasma being collisionally coupled. The fluids are described by two separate sets of non-dimensional equations 55 , which have been derived from previous models 42,56,57 . The equations are solved throughout the domain through a fourth-order  central difference scheme, and the physical variable updates  are computed by a four-step Runge-Kutta 58 scheme for time integration. The neutral fluid is governed by compressible inviscid hydrodynamics equations: while compressible inviscid resistive MHD equations model the plasma: The subscripts p and n in both plasma and neutral equations identify physical quantities of the ion-electron plasma and the neutral fluid respectively. The variables v, p, ρ, T and e are respectively the fluids velocity, gas pressure, density, temperature and internal energy, γ = 5/3 is the adiabatic index and B is the magnetic field. The terms D, R and E are respectively the source terms for mass, momentum and energy transfer between the two species, and are defined as follows: Both fluids are subject to the ideal gas law. The factor of 2 in Equation (12) is included to account for the electron pressure in the plasma contribution. The two-fluid collisional coupling is determined by the parameter α c (T n , T p , v D ), whose non-dimensional expression 59 , that includes charge exchange 60 , is found in Equation (16): where α c (0) is the initial coupling and v D = | v n -v p | is the magnitude of the drift velocity between the neutral components and the hydrogen plasma. When the magnitude of 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 expression for α c in Equation (16) takes into account the higher number of collisions occurring at supersonic drift velocities. The collisional coupling between ions and electrons is modeled by imposing a spatially uniform, constant diffusivity (η) in the system. The two fluids are in an initial thermal equilibrium, but there is not an imposed initial ionisation equilibrium. The Hall effect is not included in this study. The terms Γ rec and Γ ion are the recombination and collisional ionisation rates for a hydrogen atom. The normalised empirical forms of the rates 45,61 are: where: The rates of collisional ionisation used in this work are based on a semi-empirical model for the ionisation of hydrogen by electron impact that assumes ionisation from the ground state 45,62,63 . These rates do not include photo-ionisation or ionisation from excited states, which are known to be important for the chromosphere [64][65][66] . Two characteristic temperatures appear in Equations (17)- (20), and are based on a physical reference electron temperature T 0 in K. T e0 is the value of T 0 converted in electron volts. The initial normalisation of the system gives a bulk sound speed of unity. However, in the two-fluid model, this does not necessarily equate to a plasma temperature of unity. The ionisation and recombination rates depend on the electron temperature (which is assumed to be equal to the plasma temperature in our model), therefore a correction factor T f is needed to ensure that the desired dimensional electron temperature is being used to calculate the rates. For our normalisation, the plasma temperature is defined as: where γ(p p + p n )/(ρ p + ρ n ) is the bulk sound speed squared and is initially equal to 1. Therefore, the factor T f is applied to ensure that the electron (plasma) temperature used in the ionisation-recombination rates corresponds initially to the reference dimensional temperature T 0 . The initial ion fraction ξ p is determined from the initial choice of the temperature T 0 through its relation with ionisation and recombination rates. In a steady state, the source term for mass D in Equation (13) is zero, leading to the relation: Rearranging the terms in Equation (22) we obtain: Ionisation and recombination rates can be expressed as a function of the temperature T and the plasma density ρ p : therefore the ion fraction can be also expressed as a function of the temperature: In case of simulations not including ionisation-recombination processes, ξ p is imposed as an initial condition.
The free parameter τ IR determines the relation of the recombination time scale with the dynamic time scale of the simulation, calculated from a characteristic length and a characteristic speed of the system (see later in the section for more details). In the (PIP) code Γ ion and Γ rec are calculated in dimensional form from Equations (17) and (18), then normalised by the recombination rate in order to obtain an initial Γ rec = 1. The parameter τ IR is imposed to re-scale Γ rec , and fix an initial recombination time scale. For example, if τ IR is set equal to 10 −3 then the initial Γ rec = 10 −3 , and recombination would occur over a time scale of 10 3 . The ionisation rate Γ ion is also normalised by the same τ IR .
The terms Φ I and A heat in Equation (8) are associated to the ionisation potential and account for radiative losses. Φ I approximates the energy removed by the system through ionisation, while A heat is an arbitrary heating term included to obtain an initial equilibrium. Their non-dimensional forms are: whereΦ = 13.6/(K B γT 0 ) ensures consistency between the normalisation of ionisation potential and the equations modelling the system. Here the Boltzmann constant is K B = 8.617 · 10 −5 eV K −1 . More details on the atomic internal structure used to estimate the ionisation potential can be found in a recent paper 49 . Both sets of equations are non-dimensionalised 55 by a reference density ρ 0 and the total sound speed c s = γ(p n + p p )/(ρ n + ρ p ), initially set equal to 1. For fully ionised plasmas (MHD cases), the initial density and pressure are set constant, uniform across the domain and equal to: where ξ p is the ion fraction, equal to 1 for fully ionised plasmas. The cases that are run in a partially ionised plasma (PIP cases) present a similar normalisation as the one for a fully ionised plasma, where the bulk physical variables are set equal to the MHD values and uniform across the domain: where ξ n is the neutral fraction. This normalisation is chosen to have the physical properties dependent on characteristic length scales, which are directly comparable to the size of the plasmoids involved in the merging. The non-dimensional collision frequency can be compared to the chromospheric dimensional values by dividing it by a characteristic dimensional time scale τ col . The parameter τ col is found from the ratio of the physical plasmoid size (see Section II A for more details) and a characteristic speed in the solar chromosphere, which in this study is chosen to be the sound speed (∼ 10 km s −1 ). Therefore, the model can be easily re-scaled to examine the properties of magnetic structures of various sizes in the solar chromosphere, from a few meters to a few hundred kilometers.

A. Initial conditions
The initial setup of plasmoid coalescence is provided by a modified force-free Fadeev equilibrium 2,67 , with the magnetic field components given by Equations (33)- (35). Despite the photospheric magnetic field not being force-free at the boundary with the convective zone, its structure is rearranged before reaching the corona as the non force-free components decay due to the action of chromospheric neutrals 68 . For this reason, it is useful to impose an initial force-free field condition to the system. The classic Fadeev equilibrium magnetic field does not satisfy the condition J × B = 0 for a force-free field. Therefore the traditional Fadeev equilibrium is modified by including a component B z : In the equations above, k = π 2 and ε = 0.5, which lead to a moderately peaked current localisation at each plasmoid centre. In the limit for ε → 0 there is a less peaked current density and a weaker attraction between the plasmoids. When ε = 0, B y reduces to a current sheet characterised by the tanh profile of the well known Harris sheet 69 . The limit ε → 1 corresponds to a peaked localization and stronger attraction forces. At the upper limit (ε = 1), the current distribution becomes the delta function. The bulk (plasma + neutrals) plasma β , defined as 2(p p + p n )/B 2 , is set equal to 0.1 for all simulations. In the two-fluid cases it is possible to define a second plasma β associated to the isolated plasma 2p p /B 2 , which from the initial condition is equal to 0.002. The initial current density distribution and magnetic field lines are shown in Figure 1. Note that these initial conditions mimic our previous paper 53 .
The Fadeev equilibrium is unstable to the coalescence instability 2 . An initial velocity perturbation in both plasma and neutral species is imposed to break the initial equilibrium, given by: where v noise is a white noise component two orders of magnitude smaller than the main perturbation, included to simulate small environmental perturbations. The smaller magnitude of the white noise prevents it from dominating the motion of the two plasmoids during coalescence, while still promoting smaller scale dynamics by breaking the symmetry of the system. All simulations were performed with the same random noise seed. The sin(x) term promotes attraction of the plasmoids, while the term dependent on y localises the perturbation to a small region around the x-axis. Further details on the chosen velocity perturbation can be found in our recent paper 53 . All simulations are resolved by 4100 × 3074 grid cells, corresponding to a cell size of ∆x = 1.95 · 10 −3 and ∆y = 2.6 · 10 −3 . The resolution has been tested in order to ensure that the current sheets are resolved with our grid. The initial separation between the plasmoids, calculated from O-point to O-point, (the centre of the blue spots in Figure 1) is equal to 4L, where L is resolved by 513 grid points. The initial plasmoids width, calculated as the distance between top and bottom edges of the separatrix, is 1.66L and determined by the initial magnetic field conditions.
The dynamics of the plasmoids merging is evaluated in a squared computational domain of x = [−4, 4] and y = [−4, 4]. We use symmetric boundaries at y = −4 and y = 4: in this configuration v x and B y change sign across the boundary, while v y and B x keep the same magnitude and sign across the boundary. The boundaries at x = −4 and x = 4 are chosen to be periodic.

III. RESULTS
The coalescence instability has been well studied in MHD 2,36,70-73 . The general behaviour is that neighbouring plasmoids inside a current sheet attract each other, a second current sheet forms in between them and reconnection occurs in it, leading to the plasmoids complete merging. Extending the model to PIP (using only thermal collisions) results in a faster coalescence and the promotion of small-scale dynamics such as sub-critical secondary plasmoids formation and fractal coalescence (see our recent paper 53 ). Here the model is further extended to account for ionisation/recombination and ionisation potential effects, which provide a more realistic description of the physical processes occurring in a partially ionised plasma.
Initially, results are presented for weakly-coupled PIP cases with elastic collisions only (NIR), ionisation and recombina- a These data are the effective values of the two-fluid parameters α c and ξ p for the single-fluid cases, which are chosen as limits for the PIP simulations. b ξ p (0) and T 0 are not independent variables. The initial ion fraction is determined by setting the value for T 0 at the beginning of the calculation.
tion (IR), and ionisation, recombination and ionisation potential (IRIP), which are compared to a fiducial MHD simulation. Figure 2 shows the time evolution of the coalescence through the current density for reference cases A1 (NIR), A2 (IR) and A3 (IRIP), corresponding to the three models for PIP. From this figure, we see major differences between NIR, IR and IRIP cases in time scale, reconnection rate and dynamics of each phase of coalescence. These results are discussed in detail in Section III A. A parameter study is then performed in Section III D to investigate the effect of varying τ IR , which determines the relative importance of ionisation and recombination rates to the collision rates, and on the initial ion fraction in Section III E. The full array of simulations is shown in Table I. Simulations listed with letter M are run in fully ionised plasmas (MHD). Cases listed with letter A are the reference PIP simulations for the comparison of PIP coupling models. Cases listed with letter B compose the survey on τ IR . Cases listed with letter C compose the survey on the initial ion fraction ξ p (0).

A. Current sheet and reconnection rate
We analyse three PIP cases (listed as A1, A2 and A3 in Table I) where different coupling terms are included in each simulation: A1 corresponds to a PIP case with elastic collisions only (NIR) and its parameters are set equal to the reference PIP cases examined in our previous work 53 ; A2 includes ionisation and recombination processes (IR); A3 includes the ionisation potential (IRIP). In the cases with ionisation and recombination processes (A2 and A3) we set τ IR = 10 −3 . By choosing this value for τ IR the ratio of the collision time to recombination time is 10 −5 , which is consistent with chromospheric rates 74,75 . Based on the initial conditions, neutralion collisions occur on timescales of ∆t = (α c ρ p ) −1 = 1, ionneutral collisions on timescales of ∆t = (α c ρ n ) −1 ∼ 0.01 and recombination on time scales of ∆t = 10 3 . The coalescence instability is shown in Figure 2 at three different stages of development. In the first phase of the process the plasmoids come together, and when the oppositely directed magnetic field lines are pushed against each other a current sheet forms in the centre of the domain. In the first row of panels (a, b and c) corresponding to the end of this first phase of coalescence, the current sheet appears as the thin red vertical region at x = 0. After the current sheet formation, reconnection begins and the current sheet slowly reduces in length, together with the decrease in size of the plasmoids. The reconnected field lines form an envelope that surrounds the two merging plasmoids and constitute the external boundary of the final plasmoid resulting from the coalescence. The second row of panels (d, e and f ) shows the plasmoid merging at these later stages, while the third row (g, h and i) displays the end of coalescence and the larger merged plasmoid.
Each row of panels in Figure 2 allows to compare the three PIP cases at similar stages of coalescence. Panels (a), (b) and (c) show the current sheet structure and early small-scale dynamics at the beginning of the reconnection phase. The current sheet thickness δ is estimated by taking the full width at 1/8 of the maximum current density J z along the x-axis (y = 0): this particular ratio is chosen to be consistent with the analysis previously performed on 2.5D calculations 53 . A first difference between the three cases is observed in δ at the stage displayed by panels (a), (b) and (c). The current sheet thickess for case A1 is δ NIR = 0.018 (∼ 9 grid points), while δ IR = 0.045 (∼ 23 grid points) and δ IRIP = 0.029 (∼ 15 grid points) for cases A2 and A3, respectively. Case A1 displays the thinner current sheet and it is also the only case where the tearing instability develops. The inclusion of ionisation and recombination processes stabilises the current sheet against the tearing instability, as the plasma density increases in the reconnecting current sheet following a burst in ionisation (panels b and c of Figure 2). Varying the type of coupling between the fluids, from elastic collisions only to the inclusion of Γ ion , Γ rec and ionisation potential, the coalescence dynamics also drastically changes. Formation and expulsion of large secondary plasmoids are observed in panel (a) of Figure 2 for case A1. No plasmoids form in cases A2 and A3.
The coalescence time scale varies when varying the type of coupling between the fluids, as shown by panels (g), (h) and (i) where the final stage of the merging is displayed. Figure  3 shows the time variation of the current density at the centre of the current sheet for simulations A1, A2 and A3. These are compared to a single fluid MHD case (listed as M1 in Table  I). M1 can be treated as the limit for a completely coupled system (α c → ∞, τ IR → ∞), where density and pressure are assumed to be equal to the bulk (ion + neutral) values. The beginning of the reconnection phase is identified in all curves in Figure 3 by the first minimum in the current density. The formation of the central current sheet begins at similar times FIG. 3. Time evolution of current density J z at the centre of the current sheet for the PIP cases A1, A2 and A3 with α c = 100. The NIR case (blue), IR case (red) and IRIP (green) are compared to a MHD case (M1, black line), included as reference for the limit case of completely coupled fluids (α c → ∞). The IR and IRIP cases are run at τ IR = 10 −3 .
for all PIP simulations, suggesting that ionisation and recombination rates effects become relevant in the coalescence dynamics only once the current sheet is formed. During the initial phase of coalescence, the temperature is roughly constant and hence Γ ion and Γ rec do not vary much. During the second phase, where compressional, Ohmic and frictional heating become substantial, the temperature variation results in greatly enhanced ionisation and recombination rates (see .  Figure 4 shows the mean plasma and neutral temperature in the current sheet of our fiducial MHD and PIP simulations as they vary with time. Initially, the temperature spikes in all simulations as the current sheet collapses, due to a combination of compressional and Ohmic heating, which dominate inside the current sheet, and frictional heating. The values for T p and T n in case A1 (NIR model) are consistent with the ones found in our previous study 53 . The large increase in both plasma and neutral temperatures results from the small value of the ion plasma β and the rarefaction of the reconnection region, where the density of both fluids is very low. The neutral rarefaction observed in case A1, consistent with the results shown in Figure 4 of our previous work 53 , is the result of a strong divergence of the neutral velocity field. This divergence reaches values of 1.6 inside the forming current sheet at t = 2 and increases to be 16.7 at t = 2.4, and occurs because the neutrals are subject to two forces, generated by the drag from the plasma and by the neutral pressure respectively. As Figure 22 of our recent paper 53 shows, the neutral pressure gradient works against the drag force to hinder the neutral inflow, but it works with the drag to promote the acceleration of the material in the outflow. As result, the neutrals are able to leave the current sheet faster than they enter, resulting in the current sheet becoming rarefied.
In cases A2 (IR model) and A3 (IRIP model) both T p and T n are much lower than case A1, and there is far less temperature difference between the two species. This happens because the strong heating in the current sheet has a large effect on the ionisation-recombination rates, which in turn act on the current sheet thickness: the high ionisation rate converts neutrals to plasma which results in the current sheet thinning less, J z is smaller than NIR cases (as shown in Figure 3), and the Ohmic heating is smaller, thus leading to a cooler plasma.
In our equilibrium plasma at t = 0, the rates of both cases A2 and A3 are Γ ion = 10 −5 and Γ rec = 10 −3 , and ionisation and recombination happen on time scales of ∆t ∼ 10 5 and 10 3 respectively; in our current sheet, ionisation happen on time scales of ∆t ∼ 0.02 for case A2 (Γ ion = 47.12) and 0.05 for case A3 (Γ ion = 19.54), while recombination occurs on time scales of ∆t ∼ 8 for both cases (Γ rec = 0.13 in A2 and 0.12 in A3). This means that, for both cases A2 and A3, Γ ion varies of about 6 orders of magnitude while Γ rec varies of about 2 orders of magnitude from the initial phase to the reconnection phase. As the ionisation rate increases, so do the energy losses due to ionisation potential. The ionisation potential term acts to remove energy from the plasma, hence lowers T p . In case A3, the cooling action of the ionisation potential contributes to further decrease both T p and T n compared to case A2. Including the physics of ionisation, recombination and ionisation potential leads to a more realistic description of the system, with a temperature variation of a factor of 3 rather than a factor of 400 as observed for the NIR model.
The end of the merging is identified for each case in Figure 3 by the current density acquiring a positive value again at later times, after the development of large negative currents during the reconnection phase. The fastest coalescence occurs in case A1 (blue curve in Figure 3). The shorter time scale depends on the onset of turbulent reconnection, as the secondary plasmoids expelled from the current sheet allows a more efficient release of magnetic flux. Turbulent reconnection leads to a larger J z magnitude: the large fluctuations occurring at t > 2.8 are produced by secondary plasmoids forming and being ejected from the current sheet. In both case A2 (red) and A3 (green) reconnection occurs laminarly, similar to the MHD simulation. The coalescence time scale shortens at the inclusion of the ionisation potential. This is shown by both the comparison of the final time in panels (h) and (i) in Figure 2 and the comparison between green and red curves in Figure  3. The shortening of the coalescence time scale occurs as the cooling of the ionisation potential results in the recombination of a large portion of ions and consequently a thinner current sheet. Figure 5 shows the reconnection rate for simulations A1 (blue), A2 (red) and A3 (green) compared to the fully ionised plasma case M1 (black), where the reconnection rate as where we chose J max to be the absolute maximum value of the current density inside the current sheet, v A the initial maximum bulk Alfvén speed and B up the initial maximum value of B y in the inflow. Comparing the trend of the reconnection rate in Figure 5 with the evolution of the current density previously shown in Figure 3, it is evident that for case A1 the reconnection rate shows far less fluctuation than J z . This comes from the definition of the reconnection rate, which depends on J max : the maximum current density inside the current sheet is not necessarily occurring at its centre at all times, and this is particularly evident in case of secondary plasmoids production. While we had fixed the calculation of J z at the centre to evaluate the type of reconnection developing in the current sheet, we now want to examine the maximum reconnection rate that is achieved within the current sheet. The mean reconnection rate of case A1, where the fluids are coupled through elastic collisions and charge exchange only, is approximately M NIR = 0.12 ± 0.05: in this case the reconnection rate displays very sharp variations during the merging following the formation and ejection of the secondary plasmoids. The stabilisation of the current sheet by the action of ionisation and recombination is reflected in the lower peak value and flatter trend of the reconnection rate for the cases A2 and A3, where secondary plasmoids do not form and whose smoother fluctuations can be compared to the MHD case. The mean reconnection rates of cases A2 and A3 are also comparable to the MHD rate, M MHD = 0.047 ± 0.003. In case A2, M IR = 0.038 ± 0.006, lower than the fully ionised case. As already shown by the current density in Figures 2 and 3, coalescence in case A2 occurs over a longer time than case A1. The longer coalescence time scale is reflected in the reconnection rate, which steadily increases with time and does not show violent fluctuations, as shown by the red curve in Figure 5. The addition of the ionisation potential in case A3 speeds up the coalescence with respect to case A2, as the current sheet thins under the action of cooling and recombination, and is more compressed by the inflow. The higher mean reconnection rate of case A3, M IRIP = 0.057 ± 0.003, is consistent with the shorter coalescence time scale observed in Figure 3.
In this Section we have extensively discussed the differences between NIR, IR and IRIP coupling models for partially ionised plasmas during plasmoid coalescence. To summarise our key results, we find that ionisation and recombination have a stabilising effect on the current sheet, which is not often included in the linear theory description. Compared to fully ionised plasmas at the same bulk density, ionisation and recombination still lead to a faster reconnection in PIP, and the main effects are observed in the second phase of coalescence (reconnection phase).

B. Oscillatory behaviours of the current sheet
Several types of oscillatory motions develop during plasmoid coalescence. Waves are produced in all simulations, and are particularly evident in case A2 (IR model, intermediate initial coupling α c = 100), as shown by panels (e) and (h) of Figure 2. Smaller scale oscillations occur in the reconnection region and along the current sheet. In this Section we describe the two main oscillatory motions that develop in the current sheet.
A first type of oscillations is linked to the increase of gas pressure in the current sheet due to the plasmoids moving closer. In our previous work we observed a regular fluctuation in the O-points separation for the MHD case 53 , while such motion was suppressed in the PIP case due to the faster reconnection rate. When ionisation and recombination are included (A2 and A3 cases) we find an intermediate situation between the regular peaks as in MHD and their complete absence in NIR cases, with an initial large fluctuation shown in Figure  6 and smaller oscillations that are damped more quickly than the MHD case. The first oscillation corresponds to the large local maximum observed at the beginning of the reconnection phase for case A2 and A3, and it is a direct consequence of the processes involved in the current sheet formation. In fact, during the current sheet formation the plasma temperature increases in the thin region between the two initial plasmoids: the higher temperature leads to an initial burst in ionisation that increase the plasma pressure at the centre of the current sheet. When the plasma pressure becomes sufficiently high, it balances the pressure of the plasmoids moving closer and the plasmoids rebound off each other. At later stages, irregular oscillations of smaller magnitude can still be observed in the IR case (A2, red curve in Figure 6), while in the IRIP case (A3, green curve) the O-points distance decreases steadily. In the reconnection phase of the majority of our simulations, a second type of oscillations is observed at later times during the merging. This type of motion is displayed in Figure  7, where the divergence of the plasma velocity is displayed in the upper half of the current sheet for a reference case (B2). Oscillations start at the jet termination shock, rapidly increas- The oscillations are spatially resolved. We measure the displacement of the current sheet along x due to the oscillations by tracking the position in time of the peak in |J z |, which corresponds to the central point of the current sheet width. In case B2, half an oscillation occurs in about ∆t = 0.1, as suggested by Figure 7. At y = 0.3 the distance between the peak at t = 2.50 (first panel in Figure 7) and the peak at t = 2.60 (third panel in Figure 7) is ∆l = 7.8 · 10 −3 , which is four times bigger than the grid size ∆x = 1.95·10 −3 . Moving towards the end of the current sheet, the oscillation amplitude increases, therefore oscillations are resolved by a larger number of grid points.
In cases with very thin current sheets, such as A1 (NIR model) and the IRIP cases B1, B7, B8 and B9 in Table I, turbulent reconnection takes place before oscillations can develop. In fully ionised plasmas, where the current sheet is thicker, the oscillatory motion can be suppressed by increasing η. We display this change of regime by running an MHD case (M2) where we increase η by a factor of 2 (η = 0.003) from the value set for case M1 (see Table I). The current sheet of case M2 is subject to a steady, laminar reconnection, where no oscillations develop. The oscillatory behaviour is therefore constrained by the diffusivity, whose value determines the evolution of the current sheet dynamics into either a turbulent process (lower η) or steady reconnection (higher η).
Recent high-resolution MHD simulations of fully ionised plasmas 76 have found localised nonlinear oscillations at the edges of current sheets, studied in the framework of the dynamics of flux rope eruption in solar atmospheric plasmas. In this recent study 76 laminar reconnection occurs for a global Lundquist number S = 2.8 · 10 3 and oscillations develop when the global S = 5.5·10 3 , while at high values of S (S = 2.8·10 4 ) plasmoids form in the reconnecting current sheet 76 . The Lundquist numbers of our simulations are calculated by using the effective Alfvén speed v A,e ∼ v out , where v out is the ion outflow speed. This particular choice for approximating the Alfvén speed, consistent with our previous work 53 , accounts for the effective density given by the partial coupling between plasma and neutrals in the PIP cases. In our configuration we observe an oscillatory behaviour for Lundquist numbers between 2.6 · 10 3 and 3.3 · 10 3 in partially ionised plasma cases, while the MHD case develops this dynamics at S = 3 · 10 3 .
In order to better characterise the oscillatory motion, we look at the shift of the current sheet vertical axis position along x at a small distance from its centre. Figure 8 shows the time variation of the plasma v x at x = 0 and y = 0.2, closer to the centre of the current sheet than the region shown in Figure 7. As no grid point lays exactly at x = 0, but two sets of grid points are located at a symmetric distance from the coordinates origin, we perform a linear interpolation of the grid points located symmetrically across the y-axis. The interpolation allows to cancel out the v x,p contributions at the same magnitude but opposite in sign at the centre of the cur-  rent sheet and identify the real oscillations that displace the upper part of the structure. The beginning of the oscillatory motion is identified in this work with the position in time of the first peak having an amplitude v x,p > 0.001. We chose this threshold value as it is twice the magnitude of the initial white noise perturbation, so that random motions along the x-axis would not be mistaken with the beginning of the oscillations. In the case of the simulation B2 oscillations begin at t = 1.55.
We calculate the period of oscillations as the mean distance between the peaks of v x,p shown in Figure 8. The oscillation periods and their error, calculated as the standard deviation of the peak separation sets of measures, are reported in Table  II. Looking at the cases in the initial survey presented in Section III A, case A2 shows a longer period, P A2 = 0.39 ± 0.07, than case A3, where the period is P A3 = 0.22 ± 0.04. A period of ∼ 0.2 has been found for all the IRIP simulations run at α c = 5, as shown in Table II, which is consistent with case A3 run at α c = 100. In the MHD case, the oscillation period is P MHD ∼ 0.44, which is comparable with the IR case A2, but approximately double than all of the IRIP cases. All the simulations are highly temporally resolved, however the cadence of the time output is longer than the time step, with an output being saved approximately every 10 3 − 10 4 iterations. As our cases are analysed post process, the chosen output might set limitations on the analysis of the period. The time outputs used in this work are ∆t = 0.1 for cases M1, A2 and A3, and ∆t = 0.05 for the remaining IRIP cases. These result in having about four points to determine a period in all the cases with the exception of case A3, where only two points define the period. However, the period length was confirmed through tests run by saving smaller time outputs, which prove that the oscillatory behaviour has the same period.
It can be suggested that the longer periods identified in both MHD and IR cases are dependent on current sheet properties. The current sheets length L and thickness δ are reported in Table II at the time t peak of the first oscillation for all the cases. A correlation can be observed between oscillation period and δ : the current sheet is thicker in M1 and A2, where the period is also longer, than the IRIP cases, where the cooling from the ionisation potential results in thinner current sheets. The oscillatory dynamics does not have any influence on the reconnection rate. This can be seen through the comparison of case B2 with a second test simulation run with the same set of parameters and resolution. This second case is run in half domain in the x direction (x = [0, 4]), with a symmetric boundary located at x = 0 that prevents the onset of oscillations. Figure 9 shows the reconnection rate M of case B2 (solid black line) and the test simulation at half domain (dashed red line) where the symmetry suppresses the oscillatory motion. Both simulations display the same reconnection rate magnitude and variation across the same time interval, with minor fluctuations occurring at t > 2.
In this Section we have discussed the onset of oscillatory motions propagating in the current sheet. Previous works on fully ionised solar atmospheric plasmas 76,77 had observed the onset of localised oscillations in reconnecting current sheets. We find that the oscillations developing at the edges of current sheets are independent of ionisation and recombination processes and more generally of partial ionisation. The oscillation period shows a direct correlation with the current sheet thickness, independently of the type of plasma chosen for the simulation.

C. Secondary plasmoids
Evidence of fractal coalescence is observed in the IRIP cases B1, B7, B8 and B9, where the central current sheet is subject to the tearing instability and secondary plasmoids are produced. Figure 10 shows the interaction of two secondary plasmoids coalescing before leaving the current sheet for case B7 in Table I. Given the overall effects of ionisation and recombination on plasmoid coalescence discussed in Section III A, it is interesting to examine the secondary plasmoids and compare their properties to the ones observed in NIR simulations of our previous study 53 . Ionisation and recombination provide additional force terms that can be analysed to investigate the equilibrium of secondary plasmoids. Therefore, we look at the force balance and the magnitude of Γ ion and Γ rec across the two secondary plasmoids appearing on the left panel of Figure 10 (t = 1.15) and see if the new terms shift the force balance when compared to NIR cases 53 .
The force contributions and balance between the total pressure gradient and the Lorentz force (J × B − ∇p) are shown for the IRIP case B7 at x = 0 for t = 1.15 in the central panel of Figure 11, compared to the current density magnitude map (top panel). The yellow and black vertical dashed lines in all panels are representative of the X-point location between the two plasmoids. 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 gradient of the plasma pressure (blue curve), the magnetic pressure B 2 /2 (magenta curve) and the y-component of the magnetic tension (B · ∇) · B (red curve). This situation is similar to the cases examined in our previous study 53 : the outer structure is characterised by an almost force-free equilibrium, while −∇p p is significant around the plasmoid centre. However, a difference is observed in the distribution of the force components in the inner structure. At the plasmoid centre we observe a magnetohydrostatic equilibrium, with both magnetic pressure and magnetic tension balancing the plasma pressure gradient, while a forcefree magnetic equilibrium is sustained at the edge, where magnetic pressure and magnetic tension balance each other. Comparing the observed structure with our previous results, we see that the core region is larger in the IRIP cases than in the NIR cases. In the IRIP case the force-free equilibrium occurs in an external thin annulus, while in NIR cases force-free equilibrium nearly entirely dominates the plasmoid structure, with the exception of a very small region at the plasmoid centre. This feature is especially evident in the larger plasmoid on the right. Γ ion (red) and Γ rec (blue) are shown in the bottom panel of Figure 11. The centre of the two plasmoids, located at y = −0.204 for the smaller plasmoid and y = 0.005 for the bigger plasmoid respectively, are indicated by the two red vertical dashed lines. Recombination is observed at very small rates along the current sheet and at the centre of secondary plasmoids, while larger fluctuations are detected in the ionisation rate. Both Γ ion and Γ rec are larger in the central part of the plasmoids, where the inner structure is in magnetohydrostatic equilibrium, while the external ring characterised by a force-free magnetic equilibrium is subject to lower change rates. While Γ rec is larger at the plasmoid centre and decreases towards the ends of the inner structure, the magnitude of Γ ion in the bigger plasmoids tends to be slightly larger at the interface between inner and outer plasmoid region than in the very centre of the plasmoid.  Figure 12 shows the energy loss rate calculated as A heat − Φ I , defined by Equations (27)- (28), across the secondary plasmoids at x = 0 and t = 1.15. The cooling from the term A heat − Φ I is stronger across the secondary plasmoids than the current sheet, as shown from the troughs in the energy loss rate around the location of the plasmoids centres (identified by the red dashed lines in Figure 12). The action of the ionisation potential, which corresponds to a neat temperature decrease, can contribute to the increase in the recombination rate observed inside secondary plasmoids (see bottom panel of Figure 11).
We have seen that Γ ion , Γ rec and the ionisation potential modify the force distribution inside secondary plasmoids compared to NIR cases 53 . We want to evaluate whether the variation in these rates change the force balance inside secondary plasmoids enough to modify their structure and interaction before they leave the current sheet. From the energy loss rate we can estimate the cooling time of the plasmoid, which is calculated by dividing the total internal energy by the term A heat − Φ I . At the centre of the larger plasmoid, the cooling time is estimated to be t ∼ 0.7. The expulsion time of the larger plasmoid from the current sheet is estimated to be t ∼ 0.2, a third of the cooling time for the same plasmoid. Other smaller plasmoids coalesce or are expelled in shorter times. Therefore, we consider the secondary plasmoids to be approximately in equilibrium during their interaction inside the current sheet. Since the cooling time is larger than the expulsion time, ionisation-recombination effects do not really influence the force equilibrium of the plasmoid internal structure. Despite having a slightly different internal distribution, these secondary plasmoids act in the same way as in the NIR cases discussed in our previous study 53 .
In summary, secondary plasmoids can form in the IRIP model when turbulent reconnection develops. The force balance inside these plasmoids is slightly changed under the action of ionisation, recombination and radiative losses. However, the plasmoids equilibrium and interaction inside the current sheet are unchanged when compared to previously investigated NIR cases 53 . More details on the onset of turbulent reconnection in IRIP simulations are discussed in Section III D.

D. Survey on τ IR
The effects of the relative importance of the ionisation and recombination rates and the collision rates are investigated through a survey on the parameter τ IR , that is varied in the interval [5 · 10 −6 , 5 · 10 −1 ]. The simulations listed in Table I with IDs B1 to B6, characterised by a diffusivity η = 1.5 · 10 −3 , are compared to simulations run at a lower η = 5 · 10 −4 that are listed with IDs B7 to B9. Preliminary tests on 1D current sheets had shown that at lower collisional coupling secondary plasmoids might still form, despite the stabilisation of ionisation and recombination on reconnection. We choose to study the coalescence for α c = 5, lower than the cases presented in Section III A, as preliminary tests have identified this collisional coupling regime to better promote the onset of nonlinear dynamics in the IRIP cases. Figure 13 shows the time variation of current density J z at x = 0, y = 0 for the set of simulations at higher η (top panel) and at lower η (bottom panel). 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 initial plasmoids. The onset of reconnection begins at similar times for all the simulations in the survey, occurring later at the increasing of τ IR as ionisation and recombination processes become more important in varying the plasma composition around the X-point. For both sets of simulations, the differences between these cases in the initial phase of the coalescence are small as the initial collisional coupling is very weak, and the central current sheet forms in nearly complete absence of collisions.
In the set run at η = 1.5 · 10 −3 (cases B1 to B6 in Table I), several differences are observed in the reconnection dynamics of each simulation. Case B1, which is run at the lowest τ IR (= 5 · 10 −6 ), is the only case at higher η where secondary plasmoids formation occurs in the central current sheet, as shown by the fluctuation of the gray curve in the top panel of Figure 13. At α c = 5 the fluids are weakly coupled, with 0.05 collisions due to happen in a unit of time: this aspect allows the plasma to evolve separately with respect to the neutral fluid, and promotes secondary plasmoids formation. In case B1 the current sheet collapse (identified as the time between t = 0 and the first minimum in the current, where the current FIG. 13. Top panel: time evolution of J z at the centre of the current sheet for cases B1 to B6, run at η = 1.5 · 10 −3 . The simulations are run respectively at τ IR = 5 · 10 −6 (gray), τ IR = 5 · 10 −5 (red), τ IR = 5 · 10 −4 (green), τ IR = 5 · 10 −3 (blue), τ IR = 5 · 10 −2 (orange) and τ IR = 5 · 10 −1 (magenta). Bottom panel: time evolution of J z at the centre of the current sheet for cases B7, B8 and B9, run at η = 5 · 10 −4 . The simulations are run respectively at τ IR = 5 · 10 −5 (red), τ IR = 5 · 10 −4 (green) and τ IR = 5 · 10 −3 (blue).
sheet is compressed the most) occurs over a period ∆t ∼ 0.9: at the imposed collisional coupling, we are expecting 1 collision every ∆t ∼ 17, therefore we consider the current sheet formation to occur in an almost collisionless regime. The collisional ionisation rate of case B1, initially 5 · 10 −8 , do not significantly increase the ion fraction during the first phase of coalescence due to the lack of collisions, therefore plasmoids can form before ionisation can occur. The low ion fraction leads to a more efficient current sheet thinning, as seen in our previous paper 53 . The ejection of the first plasmoid leaves the current sheet unstable, and further smaller plasmoids are produced, leading to turbulent reconnection that efficiently reduces the coalescence time scale. Increasing τ IR by one order of magnitude (case B2, red curve in top panel of Figure 13), the ionisation rate converts a portion of neutrals large enough to increase the current sheet thickness and prevent the formation of further dynamics. For τ IR ≥ 5 · 10 −5 the tearing instability is suppressed, as the ionisation rates is sufficiently high to stabilise the current sheet.
In the simulations run at η = 1.5 · 10 −3 the magnitude of J z rapidly increases with τ IR , while the time scale for the plasmoid coalescence becomes considerably shorter. The faster coalescence is proven by the end time of the evolution in the top panel of Figure 13, where the merging completion can be identified by the current density reaching positive values (at t > 5 for all cases with the exception of case B1, corresponding to the gray curve). The reconnection phase shortens with the thinning of the current sheet. The thinning is faster for cases where Γ ion is bigger, as the cooling action of the ionisation potential is proportional to the ionisation rate. Increasing τ IR , and consequently Γ ion , reconnection becomes faster and leads to larger current densities, while the current sheet becomes thinner. In perspective, at τ IR > 0.5 we might expect a regime where the current sheet would be thin enough to promote secondary dynamics again.
In the interval of values τ IR = [5 · 10 −5 , 5 · 10 −1 ], the onset of the tearing instability might be achieved when a smaller diffusivity is chosen for the system, as thinner current sheets would form. For this reason we investigate three PIP cases (B7 to B9 in Table I) run with τ IR = [5 · 10 −5 , 5 · 10 −3 ] for a diffusivity η = 5 · 10 −4 . Bottom panel of Figure 13 shows the current density at the centre of the current sheet for the cases at τ IR = 5 · 10 −5 (red), 5 · 10 −4 (green) and 5 · 10 −3 (blue). Secondary plasmoid formation is observed for all three cases, as the smaller resistivity leads to thinner current sheets and promotes the onset of the tearing instability. At the lowest τ IR this non-linear dynamics shortens the coalescence time scale, as seen by the comparison of the red curves between top and bottom panel of Figure 13. Figure 14 shows the Lundquist number of the simulations in the survey, calculated at the beginning of the reconnection phase by using the effective Alfvén speed v A,e . In the cases at higher η, S decreases at the increase of τ IR , increasing again for τ IR > 0.005. We identify a threshold for the critical Lundquist number laying in the interval S = 5.1 · 10 3 − 7.9 · 10 4 . The simulation at the lowest τ IR (B1) is the only case at η = 1.5 · 10 −3 where secondary plasmoids are seen to form in the central current sheet. The Lundquist number for this simulation is 1.07 · 10 4 , above the critical Lundquist number. The simulations run at η = 5 · 10 −4 develop plasmoid formation at Lundquist numbers (S = 9.3 · 10 3 for B7, S = 8.2 · 10 3 for B8 and S = 7.9 · 10 3 for B9) that are consistent to the critical numbers found in our previous paper 53 .
At higher τ IR , despite having a copious production of secondary plasmoids as shown by the the fluctuation in the current density, the coalescence time scale is less affected by plasmoid dynamics. This effect might be explained by the reconnection rate saturation already observed in many 2D MHD simulations 30,42,76 , whose study revealed that above the critical Lundquist number the reconnection rate becomes almost independent of S, following the onset of nonlinear dynamics. Figure 15 shows a comparison between two PIP cases at τ IR = 5 · 10 −3 that are run at two different η, for three times (t = 1, 3 and 5). These cases, listed as B4 and B9 in Table  I, are identified by the blue curves in both top and bottom panel of Figure 13. As shown by the direct comparison of the coalescing plasmoids size between the two cases, the coalescence proceeds at a similar time scale even though in one case turbulent reconnection occurs. The Lundquist numbers of the two simulations are respectively S B4 = 2.59 · 10 3 and S B9 = 7.89 · 10 3 : secondary plasmoids in case B9 develop for an S consistent to the Lundquist numbers at which secondary plasmoids form in PIP of our previous paper 53 . The reconnection rate of both simulations is shown in Figure 16. Despite having a difference in the Lundquist number, the two rates evolves in a similar way, with the larger fluctuations occurring in case B9 following the secondary plasmoid dynamics. The mean reconnection rates, calculated between t = 1 and t = 5, are respectively M B4 = 0.062 ± 0.004 and M B9 = 0.05 ± 0.02, where the errors are calculated as the standard deviation. The two rates are close, which is consistent with the effect of reconnection rate saturation observed in other works 30,42,76 . From the Sweet-Parker steady state reconnection model, the change in η of a factor of 3 between cases B4 and B9 is expected to result in a small change of ∼ 1/ √ 3 in the reconnection rate. Re-scaling the black curve in Figure 16 by this factor to make a prediction on the reconnection rate at lower η we see a large difference of the predicted value with M B9 : the consistency between M B4 and M B9 is therefore not due to the small change in η leading to similar rates, but on the reconnection rate saturation itself. occurring once the turbulent reconnection is set.
The time variation of the mean plasma and neutral temperatures and the mean ionisation and recombination rates inside the current sheet is shown in Figure 17 for cases B4 (dashed lines) and B9 (solid lines). Γ ion and Γ rec (top panel of Figure  17) increase in both cases in the ideal phase of coalescence following the initial ionisation burst at the formation of the current sheet, and are maintained approximately constant in time in the reconnection phase, with small fluctuations around the average value. The same behaviour is observed in the temperature variation (bottom panel of Figure 17). In case B4 T p (blue dashed line) is approximately constant over the reconnection phase, while T n (red dashed line) shows small regular fluctuations at t > 2.5. Plasma is colder than the neutral counterpart. This happens because reconnection leads to a high plasma temperature initially, and once ionisation occurs large amounts of plasma energy are lost and the plasma fluid cools.
As the ionisation burst happens over a short period of time, the time scale isn't long enough for the plasma temperature to couple to the neutral temperature. In case B4, where laminar reconnection occurs, a balance is maintained with the species entering the current sheet from the inflow and the mean temperatures vary slowly, thus leading to relatively steady ionisation/recombination rates.
The fluctuations increase with the onset of the tearing instability in the current sheet as shown by the solid lines in the top panel of Figure 17. In case B9, after an increase of temperature at the current sheet initial formation, larger fluctuations appear at the onset of tearing instability. At t > 1 the plasma temperature (solid blue line) is larger than case B4, while T n is smaller. The larger T p is the result of the larger compression of the current sheet which is thinner in case B9. After t ∼ 3.5 the two mean temperatures in the current sheet reach an equilibrium and fluctuate around a similar average value. The mean temperatures at t > 3.5 are respectively T p = 1.7 ± 0.2 and T n = 1.5 ± 0.2, where the ranges are calculated as the standard deviation of the measures. At later stages there is also a general decrease of both Γ ion and Γ rec in case B9, while in case B4 the rates remain constant. The large fluctuations in presence of secondary plasmoids dynamics depend on the different mechanism of expulsion of the plasma from the unstable current sheet. The neutral temperature decreases at later times in correspondence of secondary plasmoids formationand consequently does Γ ion -and equals the plasma temperature, whose magnitude is comparable to the case without tearing instability (B4).
In this Section we investigate the relative importance of Γ ion , Γ rec to the collisional coupling α c on determining the type of reconnection during plasmoid coalescence. Because of the nature of the coalescence instability and the large compression of the current sheet in a small time interval, very small ionisation/recombination rates (τ IR ∼ 5 · 10 −5 ) are capable to affect the reconnection dynamics, as they are very sensible to temperature changes. Turbulent reconnection in weakly coupled plasmas is promoted either for nearly negligible ionisation/recombination rates (τ IR ∼ 5 · 10 −6 ) or by reducing η, which directly affects the Lundquist number.

E. Survey on the ion fraction
We investigate the changes in the coalescence dynamics of the IRIP cases following a variation of the initial ion fraction ξ p . The parameter ξ p depends directly on the reference temperature T 0 , selected at the beginning of calculation. Therefore, we evaluate changes in the coalescence instability by progressively increase T 0 by 1000 K. We compare four calculations, listed as C1, B2, C2 and C3 in Table I, where we vary the reference temperature in the range T 0 = [9855, 12855] K. Figure 18 shows the time variation of J z at the centre of the current sheet for the four IRIP cases. The changes in T 0 correspond to an ion fraction variation in the range ξ p = [2 · 10 −3 , 10 −1 ]. Note that in all cases we consider a medium that is initially dominated by neutrals.
The coalescence time scale is drastically reduced at the de-  crease of temperature, and consequently ion fraction. The time scale shortening, involving both ideal phase (initial attraction of the plasmoids to each other) and reconnection phase, is explained by the variation of the effective Alfvén speed, which scales as 1/ ξ p . The increased Alfvén speed allows more flux to enter the reconnection region, hence feeding the reconnection process and accelerating it.
In terms of the observed trend in the development of the current density magnitude, this survey reproduces the results of our previous paper's study on ξ p 53 . Going towards higher ξ p , where the magnetic forces are felt by larger portions of the fluid, the coalescence time scale tends to the MHD case. The reconnection rate increases slightly as ξ p (0) increases. This trend, shown by the values of mean reconnection rate collected in Table III, is opposite of what was observed in our previous paper 53 . The reversal in trend might be explained by examining both the variation in plasma and neutral densities across the current sheet, and the average ionisation and recombination rates. Figure 19 shows the mean values of Γ ion and Γ rec (top panel) and of ρ p and ρ n (bottom panel) across the current sheet for the four cases presented in the survey. The oscillations of the average value of ρ p and ρ n can be associated to the oscillatory motion of the current sheet.
The compression of the current sheet by the merging plas-  18. Time evolution of current density J z at the centre of the current sheet for four PIP cases (C1,C2,C3 and B2). The initial reference temperature T 0 is respectively T 0 = 9855 K, corresponding to an initial ξ p = 2 · 10 −3 (turquoise), 10855 K, corresponding to ξ p = 10 −2 (red), 11855 K, corresponding to ξ p = 4 · 10 −2 (purple) and 12855 K, corresponding to ξ p = 10 −1 (orange). The current density of an MHD case (M1) is included as comparison to the limit T 0 → ∞. moids results in heating that leads to a burst in ionisation, observed in particular at the formation of the current sheet. At the increase of the initial ξ p , both Γ ion and Γ rec decrease inside the current sheet, as shown in the top panel of Figure 19. This is due to the larger ρ p at the centre of the structure, as the plasma pressure prevents the current sheet to thin more, hence showing a smaller temperature change. When the ionisation rate is small it leads to a limited ionisation, which also turns into a lesser thinning of the current sheet as the cooling action of the ionisation potential is reduced. The lower cooling affects the recombination rate, whose decrease leads to fewer  Having an initial lower ξ p , and consequently a larger ionisation rate due to the initial current sheet thinning, both ρ p and ρ n increase more inside the current sheet, the latter due to the larger cooling factor provided by the ionisation potential. Although the coalescence is faster for lower ξ p (t = 0), the local increase in density leads to a comparable reconnection rate for all the cases in this survey.
The change of the effective Alfvén speed and plasma density affects the Lundquist number, which in turn extends the time scale of the reconnection phase. The values for the Lundquist number calculated at the beginning of the reconnection phase are shown in Table III, and it is shown that the Lundquist number also slightly increases in the interval. The value of S for case C3 (ξ p = 10 −1 ) is consistent with our previous study without ionisation and recombination 53 . At lower ξ p (t = 0) the Lundquist number is generally lower than our previous work, but this is due to the increase in plasma density that comes from the ionisation at the centre of the current sheet.
In general, the evolution of the coalescence instability observed from changes in the current density follows the same trend as the NIR cases previously studied 53 . However, ionisation/recombination rates act on both the reconnection rate M and the Lundquist number S. Both parameters vary less than the NIR simulations in the same ξ p (0) interval.

IV. DISCUSSION
Plasmoid coalescence is a very important process for promoting fast magnetic reconnection in many reconnecting systems: in fact, plasmoids allow the shortening of the time scale for explosive events by acting on the size of reconnecting current sheets. It is not entirely clear about how coalescence develops in a partially ionised plasma and in particular how it is affected by the charge-neutral species interaction. In this work we have extended our previous study on the coalescence instability in partially ionised plasmas 53 by including ionisation, recombination and optically thin radiative losses in the two-fluid coupling terms. The goal of our update is to model a more realistic coupling between different particle species during coalescence. We have observed several changes in the development of this instability that are related to the newly included two-fluid physics: 1. Ionisation and recombination have a stabilising effect on the current sheet as they are responsible for thickening it, but still lead to a faster reconnection in PIP when compared to MHD cases at the same bulk density.
2. The formation of the central current sheet begins at similar times for all PIP cases at the same initial collisional coupling, with or without ionisation and recombination, ionisation potential and background heating. This is because in the first phase of coalescence the rates are in equilibrium, and the temperature on which they depend does not show sharp variation. During the reconnection phases, where the temperature increases, ionisation/recombination effects dominate the reconnection dynamics, and are responsible for the type of reconnection that develops in the current sheet.
3. The internal structure of secondary plasmoids is slightly altered under the action of ionisation and recombination: a larger part of their structure is in magnetohydrostatic equilibrium, while the external region characterised by a force-free equilibrium is reduced to a thin annulus. However, the secondary plasmoid dynamics inside the current sheet is unchanged when compared to NIR cases 53 .
4. Small initial ionisation/recombination rates (∼ five orders of magnitude smaller than the collision rates) in weakly ionised plasmas can still lead to an efficient stabilisation of the current sheet against the tearing instability in the reconnection phase. This happens because the temperature increases locally enhance the rates by several orders of magnitude inside the current sheet, thus leading to a substantial increase in ρ p inside the current sheet.
5. When varying the initial ion fraction, the general development of the coalescence is unchanged from NIR cases, but ionisation and recombination influence reconnection rate and Lundquist number because they increase locally the plasma density, thus affecting the value of the effective Alfvén speed.
In addition to the results that are directly dependent on the presence of ionisation/recombination rates, we also find that in cases of laminar reconnection the current sheet develops an oscillatory behaviour in the outflow that is independent of partial ionisation effects, and correlated to the current sheet thickness (see Figure 7). From our findings, the stabilisation factor in (1) provides interesting consequences for a recent paper 44 that had investigated fractal reconnection in partially ionised plasmas in strongly, intermediate and weakly coupled regimes. As changes in the ionisation level directly affects the current sheet evolution, then the time scales for fractal tearing previously identified 44 will be largely affected by the presence of ionisation-recombination processes. The local dynamics observed in (1) and (2) is also consistent with the ionisation bursts found in previous studies on magnetic reconnection in multi-fluid partially ionised plasmas 7,42,43,[46][47][48] .
We conclude with the implications of our findings for chromospheric reconnection. The chromospheric rates of ionisation and recombination are typically 74,75 in the range 10 −5 −10 −3 s −1 , between 2 and 7 orders of magnitude smaller than the typical chromospheric neutral-ion collision frequency (10 −1 − 10 2 s −1 ) 78-80 . Our study, where rates are set consistent with the chromospheric values, demonstrates that ionisation and recombination rates are large enough to suppress the onset of fractal coalescence and small-scale dynamics, as they act on current sheets properties. However, multi-fluid physics is still capable to promote fast reconnection, hence explaining the short time scales of chromospheric explosive events. A process that is not accounted for in this work is the interaction of the plasma with the photospheric radiation field. While collisional ionisation by electrons is important for heavier elements, photo-ionisation is the predominant process for ionising hydrogen in the chromosphere 64 and leads to a hydrogen ionisation rate 65,66 Γ ph ion ∼ 1.4 · 10 −2 s −1 that can become orders of magnitude bigger than collisional ionisation by electrons 66,81 Γ ion ∼ 7.8 · 10 −5 s −1 . The inclusion of photo-ionisation would impose changes in the balance between ions and neutrals as function of the atmospheric temperature, therefore it will be subject for further updates of our model.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request. The (PIP) code is available at the following url: https: //github.com/AstroSnow/PIP. Details of the code and equations are available in A. Hillier, S. Takasao & N. Nakamura, A&A 591, A112 (2016) 54 . A copy of the code version used in this work and the initial condition of the reference simulations can be found in a dedicated repository provided by G. Murtas: https://github.com/GiuliaMurtas/PIP.git.